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I.  Technical  Report  Summary 

Two  major  advances  have  been  made  in  the  research  on  synthesis  of  long- 
period  S  waves  at  distances  less  than  40°:  the  incorporation  of  inelastic 
attenuation  and  the  development  of  a  code  to  include  a  representation  of 
the  source  function.  The  large  effects  on  the  waveforms  of  rapid  variations  of 
velocity  in  the  upper  mantle  can  be  calculated  and  separated  from  the 
effects  of  anelasticity . 

Examples  of  synthetic  seismograms  are  included  to  illustrate  the  effects 
of  attenuation  on  the  displacement  waveforms,  as  well  as  the  effects  of  the 
response  of  a  typical  long-period  seismograph  system. 

The  application  of  Green's  Function  techniques  to  elastodynamics  has 
led  to  methods  for  treating  a  variety  of  problems  in  wave  propagation  aid 
earthquake  source  representations.  The  complete  theory  is  presented  in 
this  report.  Wave  propagation  is  a  realistic,  layered  earth,  generalized 
by  a  nonlinear  source  and  the  dynamic  field  due  to  stress  relaxation  around 
a  geometrically  general,  growing  inclusion  (an  earthquake  source)  in  a 
spatially  heterogeneous  initial  stress  field  are  two  of  the  significant 
problems  that  have  been  treated.  Computations  based  on  the  equations 
derived  here  are  currently  being  carried  out. 

There  is  evidence  that  specific  anelastic  attenuation,  expressed  as  Q,  is 
frequency  dependent.  Q  models  of  the  Earth  based  on  free  oscillations, 
surface  waves  and  body  waves  show  that  in  the  broad  frequency  band  covered 
by  the  input  data,  Q  increases  with  frequency.  Frequency-dependent  Q  is 
modelled  by  a  relaxation  process  with  a  range  of  relaxation  times.  An 
investigation  of  the  relaxation  time  characterizing  the  high  frequency 
corner  of  the  absorption  based  was  carried  out  using  the  data  from  21 
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shallow  earthquakes,  1  at  intermediate  depth,  and  4  deep  ones. 

The  criteria  used  was  that  the  Q-corrected  spectrum  show  decay  at  high 

-2  -3 

frequencies  at  a  rate  in  the  range  f  to  f  .  The  depth  dependence 
of  the  resulting  relaxation  times  and  corresponding  values  of  T/Q 
was  examined.  A  mixed  effect  of  depth  and  frequency  dependence  of  P-wave 
attenuation  was  found.  The  important  conclusions  is  that  the  P  and  S 
attenuation  data  can  only  be  reconciled  by  including  a  bulk  loss  mechanism, 
in  addition  to  a  shear  loss  mechanism.  Although  the  results  are  not 
unique,  this  suggests  that  the  bulk  loss  mechanism  is  operative  in  the 
upper  mantle,  perhaps  within  the  asthenosphere. 


II.  Synthesis  of  Long  Period  Body  Waves  in  an  Anelastic  Earth 

V.F.  Cornier 

Sunthesis  of  long  period  SH  body  waves  has  been  completed  in  the  PEM-C 
(Figure  2.1)  earth  model  of  Dziewonski  et  al.  (1975).  This  earth  model 
in  conjunction  with  a  simple  frequency  independent  Q  model  was  chosen  as  a 
starting  model  for  inversion  for  the  upper  mantle  structure  of  western 
North  America.  Synthetic  modeling  of  waveforms  at  distances  less  than  40 
degrees  must  be  undertaken  to  separate  the  large  effects  of  rapid  velocity 
variation  in  the  upper  mantle  from  the  effects  of  anelasticity .  The  in¬ 
corporation  of  attenuation  in  the  analysis  through  a  frequency  dependent 
complex  velocity  profile  has  been  taken  as  complementary  check  of  spectral 
studies  of  upper  mantle  attenuation  structure. 

Seismograms  have  been  synthesized  using  the  method  described  by 
Cormier  and  Richards  (1977),  in  which  the  Fourier-transformed  displacement 
is  evaluated  by  the  evaluation  of  an  integral  over  short  paths  in  the 
complex  ray  parameter  plane.  Langer's  approximation  to  the  radial  eigen¬ 
functions  in  the  integrand  of  reflection-transmission  coefficients  corrects 
for  the  effect  of  velocity  gradients  and  boundary  curvature  of  inhomogeneous 
spherical  layers.  In  the  last  funding  period  a  significant  improvement  in 
the  computation  speed  of  synthetics  has  been  achieved  by  using  the  three 
point  integration  formula  described  by  Jeffreys  and  Jeffreys  (1956)  for 
the  evaluation  of  the  quantity 


in  inhomogeneous  layers.  This  quantity  is  required  to  evaluate  the  Langer 
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approximation  (Richards,  1976).  By  use  of  the  three  point  formula,  the 
computation  of  synthetic  seismograms  for  direct  SH  body  waves  in  an  attenuative- 
dispersive  earth  model  with  two  upper  mantle  T-A  triplications  requires 
only  10  minutes  execution  time  on  a  CDC-6600  at  2  degree  intervals  from 
20  to  34  degrees. 

Multiple  reflections  along  the  underside  of  discontinuous  velocity 
increases  in  the  upper-mantle  of  PEM-C  are  included  by  use  of  generalized 
reflection  coefficients  and  integration  paths  of  the  type  described  in 
Cormier  and  Richards  (1977).  For  example,  the  reflection  coefficient  at 
critical  incidence  is  taken  as  the  functions 


.(2) 


22  (1)  » 
°2 


,(2) 

2 

,C1) 


„(2) 

0  2  °1 


21  12  (2)  (2) 
°2  °1 


.(1) 


1  -  R 


11  (2) 

°1 


along  the  integration  paths  illustrated  in  Figure  2.2,  where  R22  is  the 
reflection  coefficient  of  SH  from  the  top  of  the  boundary,  R^  the 
reflection  coefficient  from  the  bottom,  transmission  coefficients 

from  the  top  or  bottom  respectively,  »  °2^  radial  eigenfunctions 

for  (1)  up  or  (2)  downgoing  SH  waves  in  the  upper  medium, 
radial  eigenfunctions  in  the  lower  medium.  These  coefficients  are  related 
by  the  equation 
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At  distances  much  greater  than  critical  a  reflection  coefficient  is  constructed 
from  the  functions  above  for  use  along  the  integration  path  illustrated 
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Integration  paths  near  a  cusp  at  critical  incidence.  The 
coefficients  are  taken  along. the. segments  as  shown.  (Ratios 
of  radial  eigenfunctions  a  fa  are  omitted  in  the  figure 
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The  N  waves  constituting  the  direct  transmitted  and  N-l  multiples 

are  included  in  a  separate  ray  parameter  integral  when  they  are  well  separated 

in  arrival  time  from  multiples  described  by  eq.  (2.2). 

By  proper  choice  of  combinations  of  these  functions  and  integration 
paths,  the  effect  of  all  the  interfering  multiples  in  two  upper  mantle 
triplications  can  be  included  in  the  evaluation  of  a  single  integral  over 
ray  parameter  in  the  distance  range  near  20°. 

Figure  2.4a-b  shows  the  results  of  the  synthesis  for  the  response 
of  the  PEM-C  model  to  a  SH  source  for  the  direct  body  wave.  Comparison 
of  the  results  with  (2.4b)  and  without  (2.4a)  attenuation  demonstrate 
that  the  removal  of  high  frequencies  by  attenuation  obscures  details 
in  the  waveform  due  to  multiple  arrivals  from  first  order  discontinuities 
in  the  earth  model.  The  difference  in  model  response  remains  visible  even 
after  inclusion  of  the  transfer  function  of  a  long  period  seismograph 

(Figures  2.5a-  .5b).  These  results  are  consistent  with  those  reported  by 

$ 

Rennet  (1975)  for  inclusion  of  attenuation  in  the  reflectivity  method  of 
seismogram  synthesis. 
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Figure  2.3:  Integration  path  for  an  interference  head  wave. 
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Figure  2.6  compares  synthetics  generated  by  convolving  the  model 
response  with  a  source  function  determined  for  the  Borrego  earthquake 
(Helmberger  and  Engen,  1974)  with  observed  seismograms  of  the  earthquake. 

The  source  function  includes  near  source  multiples  and  attenuation  through  the 
Carpenter  (1966)  Q  operator.  The  model  response  for  no  attenuation 
consequently  was  the  one  convolved.  Future  research  will  compare  wave¬ 
forms  predicted  with  the  Carpenter— Q  operator  with  those  determined  by  _  ; 
including  attenuation  in  the  frequency  domain.  The  latter  approach  has  the 
advantage  of  more  easily  handling  complicated  frequency  dependence  of  Q, 
but  requires  slightly  longer  computation  time. 

Analysis  of  Figure  2.6  indicates  that  of  the  later  arrivals  predicted 
by  the  PEM-C  model,  the  one  due  to  reflection  from  the  discontinuity  at 
671  km  is  not  greatly  in  disagreement  with  the  observed  seismograms.  The  arrival 
forming  the  latest  and  largest  trough  in  the  synthetics  corresponds  to  a  wave 
diffracted  along  the  top  of  the  lower  boundary  of  the  low  velocity  zone  of 
PEM-C.  Revision  of  the  PEM-C  model  to  agree  with  waveform  (Figure  2.6) 
travel  times  (Figure  2.7)  can  be  simply  accomplished 

by  modifying  the  structure  of  the  low  velocity  zone.  Unfortunately  wave¬ 
forms  cannot  sufficiently  constrain  deeper  structure  because  a  long  period 
recording  of  a  sufficiently  large  earthquake  cannot  resolve  the  features 
of  waveform  due  to  interfering  triplications  near  20  degrees.  Attenuation 
makes  this  resolution  more  difficult.  One  must  rely  on  carefully  matching 

amplitudes  near  20°  and  matching  high  quality  waveforms  at  distances  greater 

$ 

than  30  degrees  uncontaminated  by  PL  waves. 

In  the  next  research  period  SH  synthesis  will  be  completed  in  earth 
models  satisfying  travel-time,  amplitude,  and  waveform  data  for  western 


URCE 


jure  2.6:  Synthetics  (30,  32,  340)  generated  by  convolving  model  response  with  Burrego  Mountain  earthquake 
(April  9,  1966)  that  given  by  Helmberger  and  Engen  (1966);  data  at  WWSSN  stations  SCP  and  OGD. 


Figure  2.7:  Reduced  travel  time  curve  in  PEM-C  compared  with  observed  travel 


North  America.  A  range  of  upper  mantle  models  for  S  velocity  and  Q  will  be 
suggested.  The  validity  of  including  attenuation  in  the  time  domain  via 
the  Carpenter  Q  operator  will  be  investigated. 
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III.  Anelastic  Properties  of  the  Earth  Inferred  from  Earthquake  Spectra: 
Investigations  of  Frequency  Dependent  Q  Models 

G.M.  Lundquist 

Introduction 

Anelastic  attenuation  of  seismic  waves  is  parameterized  by  the  specific 
quality  factor,  Q  ,  defined  by 


_1_  AE 
2  it  E 


(3.1) 


where  AE  is  the  energy  dissipated  in  a  cycle  whose  average  energy  is 
E  .  In  a  layered  medium,  the  wave  is  attenuated  proportional  to  a  weighted 
average  of  the  Q  ^  in  each  layer,  where  the  weight  is  given  by  the  travel 
time. 


(3.2) 


where  t^,  is  the  travel  time  in  each  layer.  The  functional  dependence  of 
the  absorption  becomes  apparent  in  the  definition, 

_  uK 

TOTAL  ANELASTIC  ATTENUATION  =  e  2QEff  (3.3) 

where  T  is  the  total  travel  time.  For  convenience,  we  will  note  the 

ratio  T/Q  ,  by  t  with  subscript  a  or  0  for  P  or  S  waves  respectively. 

i 

* 

The  minimum  functional  dependence  for  t  must  be  on  velocity  and  Q  structure 
along  the  ray  path. 

If  an  additional  dependence  upon  frequency  is  to  be  found,  then  the 
velocity  and  Q  structures  must  either  be  known  or  somehow  eliminated  from 
the  problem.  As  before  (Lundquist,  1975)  we  assume  a  multiplicative 
separability  between  the  depth  dependence  and  frequency  dependence  of  the 


form 
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t  (r,  Cd)  =  t*(r)  R(u>)  (3.4) 

where  r  is  radius  in  a  spherically  symmetric  earth  model  and  R(oj)  carries 

all  of  the  frequency  dependence.  Equation  (3.4)  demonstrates  explicitly 

the  nonuniqueness  inherent  in  determining  R((jj)  when  t  is  not  known. 

* 

Unfortunately,  since  t  appears  in  an  exponent,  ratioing  techniques  will 
not  isolate  R(w)  .  That  is,  the  base  t  functions  must  be  estimated 
independently. 

Toward  this  end,  a  review  has  been  made  of  published  velocity  and  Q 
models.  Rather  than  finding  a  consensus  among  Q  models,  a  frequency 
dependence  is  noted  depending  upon  the  data  set  used.  In  general,  Q  seems 
to  be  an  increasing  function  of  the  frequency  of  the  data  set.  To  determine 
whether  this  frequency  dependence  may  be  related  to  R(u>)  as  determined 
from  body-wave  spectra,  preliminary  attempts  are  made  to  model  the  differences 
between  published  Q  structure  as  a  function  of  a  frequency  dependent  Q, 
with  quite  good  results. 

Published  Models. 

Velocity  models  are  quite  well  constrained  by  both  body-wave  and  free- 
oscillation  inversions.  The  slight  (1%)  base-line  shift  in  theoretical 
travel  times  was  shown  to  be  a  result  of  ignoring  anelasticity  in  the  free 
oscillation  studies.  Hart  et  al.  (1977)  corrected  the  entire  set  of  known 
spheroidal  and  toroidal  modes  for  attenuation  and  reduced  the  base-line 
shift  from  the  Jeff reys-Bullen  tables  to  less  than  a  second.  Since  the 
resulting  velocity  model,  QM2,  agrees  with  both  body-wave  and  free-oscillation 
data,  it  has  been  adopted  for  use  in  the  present  study. 

Q  models,  on  the  other  hand,  are  neither  as  well  constrained  nor  as 
consistent.  The  range  of  models  is  shown  in  Figure  3.1.  Model  SL1  was 


DEPTH (km) 


Figure  3.1:  Q  vs.  depth.  Models  SLl  and  MM8  were  derived  from  free  oscillations 
and  surface  waves,  respectively.  Model  AFL  was  derived  from  body 


waves. 
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derived  from  free-oscillation  data  (Anderson  and  Hart,  1977);  MM8  was 
derived  from  surface-waves  (Anderson  and  Archambeau,  1964);  and  AFL  was 
derived  from  P  body  waves  (Archambeau  et  al.,  1969).  The  frequencies  in 
the  data  sets  appropriate  to  each  model  are  given  in  Table  3.1.  Note 
that  while  the  free-oscillation  and  surface-wave  models  overlap  in  frequency, 
the  body-wave  model  is  taken  from  waves  in  a  completely  independent  frequency 
range . 

TABLE  3.1 


Model 

Data  Frequencies 

SL1 

.0003-. 015  Hz 

MM8 

.0025-. 02  Hz 

AFL 

2-5  Hz 

The  variation  in  Q  models  was  examined  as  a  function  of  t 

* 

Figure  3.2  shows  t  vs  epicentral  distance  for  a  shallow  focus  earthquake 

for  each  of  the  models  of  Figure  3.1.  Model  SL1  included  both  Qa  and 

Qg  .  For  the  other  two  models,  Qg  was  obtained  by  =  2.35  Qg 

When  combined  with  velocity  model  QM2,  this  Q  ratio  gives  ta/tg  =  4.3  , 

*  * 

while  the  ratio  for  SL1  is  t  /tQ  -  4.55  .  There  is  very  little  change  in 

CL  p 

* 

the  shape  of  t  vs  distance  as  a  function  of  source  depth  for  any  particular 
model . 

In  both  Figures  3.1  and  3.2,  the  low-frequency  models  overlap,  as  might 

* 

be  expected  since  they  deal  with  data  in  the  same  frequency  range,  t  (AFL)  , 

* 

however,  is  only  half  of  the  values  for  t  (MM8)  because  Q(AFL)  is  con¬ 
sistently  higher  than  that  for  the  other  models.  If  both  SL1  and  AFL  are 
correct,  then  Q  must  increase  with  frequency. 


T/Q  (Sec) 


t*  vs.  DISTANCE 
VELOCITY 
MODEL  QM2 


Qig  from  SLI 


O^g  from  MM8 


Qp  from  AFL 


JO  20  30  40  50  60  70  80  90 

EPICENTRAL  DISTANCE  (Deg) 

Figure  3.2:  T/Q  vs.  Distance.  The  velocity  model  used  was  QM2.  Note  that 
T/Q(AFL)  is  well  below  T/Q  for  the  low  frequency  models. 
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Figure  3.3  shows  the  variation  of  t  vs  hypocentral  depth  as  a 
function  of  percent  of  surface  focus  value.  Low  frequency  waves  apparently 
miss  about  23%  of  the  total  attenuation  if  the  source  is  at  600  km  depth. 

High  frequency  waves,  on  the  other  hand,  apparently  miss  more  than  40%  for 
the  same  change  in  source  depth,  and  the  difference  seems  to  occur  over 
the  depth  range  of  the  asthenosphere .  This  difference  may  also  be  reconciled 
simply  as  a  function  of  frequency  if  low  and  high  frequency  waves  see 
different  changes  in  apparent  Q  as  a  function  of  source  depth. 

Q(f )  Model. 

Before  going  further,  it  is  appropriate  to  briefly  review  the  form  of 
the  frequency  dependence  of  Q  used  here.  The  function,  R(w)  ,  in  the 
attenuation  exponent  of  equation  (3.3)  generates  an  absorption  band  from 
a  constant  T/Q^^.  •  The  absorption  band  is  constructed  theoretically 

as  the  superposition  of  specific  absorption  mechanisms  each  of  which  is 
modelled  by  a  standard  linear  elastic  solid.  Each  separate  mechanism 
attenuates  according  to 


Q  1  =  C 


WT 


i  ,  2  2 

1  +  U)  T 


where  C  is  a  constant  depending  upon  the  elastic  parameters  of  the 

medium  and  T  is  the  medium  relaxation  time  (see  Mascn,  1958).  A  distribution 

of  relaxation  times  of  the  form 


D(r) 


(  1/x  ,  >  T  >  X2 

*  °  ,  T  >_  T1  T  <_  T2 


defines  an  absorption  band  as  (Liu  et  al.,  1976) 


ras%  SURFACE  FOCUS 
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Figui a  3.3:  T/Q  vs.  Source  Depth.  The  velocity  model  used  was  QM2.  Note 
the  variation  as  a  function  of  the  frequency  in  the  data  set 
used  to  generate  the  Q  model. 
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The  parameter  of  interest  in  (3.5)  is  x^  ,  the  high-frequency  half¬ 
power  point  of  the  absorption  band.  The  other  corner  period  is  arbitrarily 
placed  well  beyond  the  longest  periods  of  interest  in  body-wave  studies  at 
=  2000  sec  .  Thus  the  only  frequency  dependence  introduced  in  this  work 
is  a  roll-off  in  absorption  toward  high  frequencies  which  may  be  varied  by 
changing  T 2  . 

The  point  of  introducing  a  frequency  dependent  Q  is  to  recover 

observed  seismic  source  spectra  from  seismograms.  Theoretical  source 

spectra  all  show  amplitude  decays  toward  high  frequencies  of  the  order 
-2  -3 

u>  to  ii)  .  The  slope  of  an  observed  spectrum,  however,  is  controlled 
much  more  by  the  anelastic  attenuation  along  the  travel  path  than  by  initial 
source  properties.  If  a  Q  correction  function  is  appropriate,  it  should 
give  back  an  observed  spectrum  with  slopes  in  the  theoretical  range. 

Thus  R(w)  is  used  in  the  attenuation  correction  as  a  decay  slope 
modification.  The  results  of  that  modification  will  be  presented  in  the 
next  section. 

Body  Wave  Spectra. 

To  date,  21  shallow  earthquakes,  A  deep  and  1  intermediate  depth  event  have 
been  studied.  Body-wave  spectra  were  computed  from  digitized  seismograms 
using  an  autocorrelation  technique  which  smoothed  the  individual  spectra. 
Additional  smoothing  was  obtained  by  averaging  several  stations  for  each 
wave  type  of  each  event  studied.  The  result  was  usually  (but  not  always) 
a  noise-free  spectrum  whose  high  frequency  decay  slope  could  be  uniquely 
fit  by  a  single  line  segment. 
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An  average  spectrum  without  an  anelastic  attenuation  correction  was 

prepared  for  each  event,  under  the  assumption  that  T/Qe^  is  nearly 

constant  with  epicentral  distance.  Then  each  average  spectrum  was  corrected 

* 

for  attenuation  with  a  range  of  •  Initially,  the  base  t  function 

*  *  * 

was  assumed  to  be  t  =1.0  and  t„  =  4.0  ;  and  the  variation  of  t 

a  p 

with  hypocentral  depth  was  chosen  to  be  that  of  model  AFL.  After  the 

* 

review  of  published  models  was  begun,  this  choice  of  t  was  seen  to  mix 

* 

high  and  low  frequency  models.  Thus  the  data  set  was  reprocessed  with  t 

~  *  * 

essentially  following  model  SL1.  For  that  test,  t^  =  1.0  and  tg  =  4.5 
for  shallow  events,  and  the  SL1  depth  dependence  from  Figure  3.3  was  applied 
for  deep  and  intermediate  events. 

The  results  of  the  tests  are  given  in  Tables  3.2  and  3.3  for  shallow 

-3 

and  deep  events,  respectively.  The  values  ~  00  for  to  slopes  of 

deep  P  waves  result  from  the  fact  that  the  uncorrected  station  spectra 

-3 

already  have  slopes  less  than  <o  .  Since  the  Q  correction  can  only 

_3 

decrease  the  decay  rate,  an  <u  decay  rate  cannot  be  obtained  for  any 

•  This  may  imply  a  significant  difference  between  deep  and  shallow 

-3 

source  mechanisms,  or  it  may  imply  that  the  assumption  of  to  slope  is 
incorrect. 

In  Tables  3.2  and  3.3,  the  values  of  in  parentheses  are  the 

* 

results  of  reprocessing  with  t  defined  by  model  SL1.  Note  that  the 

* 


result  of  increasing  t  is  to  increase  the  applied  attenuation  correction 

and  raise  the  slope.  Thus  a  larger  is  required  to  get  back  the  same 

•k 

slope  for  increased  t  .  Though  the  difference  is  negligible  for  shallow 

* 

S  waves,  the  difference  between  t  from  models  AFL  and  SL1  (or  MM8) 

as  a  function  of  source  depth  is  significant.  If  any  physical  understanding 


i 
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Date 

72/08/30 

72/08/30 

72/08/09 

72/04/09 

71/05/22 

71/04/03 

71/03/23 

70/08/13 

70/06/05 

70/03/28 

69/02/11 

67/08/15 

67/02/11 

63/04/19 

71/01/10 

70/07/26 

69/11/07 

69/02/28 

65/06/27 

65/02/02 

64/03/28 

T/Qeff(P) 

T/Qeff(s) 

T/Qeff(s) 


TABLE  3.2 

EVENT  LIST  SHALLOW  EARTHQUAKES 


Depth 

T,  (w-2) 

2  P 

T,  (<xf3) 

2  P 

T2  s<“"2) 

T2 

33 

.09 

.17 

.17(.18) 

,24(.24) 

33 

.11 

.18 

.16(. 16) 

,24(.24) 

15 

.06 

.12 

33 

.08 

.16 

.17 ( .18) 

.26(. 26) 

33 

.08 

.18 

. 17(. 18) 

,26(.26) 

33 

.05 

.13 

.19(.20) 

.28 (.28) 

33 

.06 

.14 

(.19) 

(.28) 

15 

.11 

.20 

20 

.08 

.18 

.18 (.18) 

.27(.27) 

15 

.11 

.20 

33 

.08 

.16 

,19(. 20) 

,29(.29) 

33 

.08 

.18 

,19(.20) 

.28(.28) 

5 

.08 

.18 

33 

.08 

.16 

,18(.18) 

.26(.26) 

33 

.05 

.14 

35 

.05 

.16 

35 

.04 

.14 

22 

.02 

.09 

27 

.09 

.20 

12 

.06 

.14 

21 

.08 

.20 

1.0 

4.0  for  values  not  in  parentheses 
4.5  for  values  in  parentheses 


TABLE  3.4 

STATISTICS  OF  OBSERVED  RELAXATION  PARAMETERS 
(t  from  Model  SL1) 

Average 


t2  p(“‘2) 

T2  p(“‘3) 

.  -2 
T2  s(b) 

Shallow  Events 

.073 

.162 

.185 

Deep  and  Intermediate 

* 

Events 

.094 

.22 

.146 

From  event  of  72/03/30  only 
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is  to  be  obtained  from  the  x^  determined  from  this  study,  the  depth 

* 

dependence  of  the  base  t  function  must  be  known. 

Added  emphasis  to  this  last  conclusion  is  given  by  Figure  3.4,  where 

* 

is  plotted  as  a  function  of  the  base  t  value.  The  purpose  of  this 

k 

test  is  multifold:  (1)  Examination  of  the  values  of  t  and  x  in 

* 

Tables  3.2  and  3.3  showc  a  montonic  trend  in  t  vs  t ^  from  deep  P  waves 
to  shallow  S  waves,  and  it  is  necessary  to  show  that  the  results  presented 
here  are  not  artifacts  of  the  processing.  (2)  The  effect  of  base  Q  model  on 
must  be  investigated.  (3)  The  deep  event  of  72/03/30  was  chosen  for 

_3 

the  study  because  it  is  the  only  deep  event  for  which  w  slopes  may  be 
obtained  in  the  P-wave  spectrum.  Since  this  event  yields  the  only  data  in  the 

_3  * 

u>  column,  it  should  be  carefully  checked.  Figure  3.4  shows  that  t  vs 

is  indeed  monotonic  increasing,  but  the  curves  are  distinctly  different 

for  P  and  S  waves.  The  change  in  T 2  as  a  function  of  wave  type  is  not  an 

* 

artifact.  There  is,  however,  an  interdependence  between  t  and  such 

that  the  change  in  T 2  as  a  function  of  source  depth  may  be  strictly  a 

* 

matter  of  the  different  t  required. 

Table  3.4  shows  that  T2(P)  increases  with  source  depth,  while 

TjCs)  decreases.  The  decrease  in  T2(S)  is  plotted  in  Figure  3.4,  and 

* 

closely  resembles  the  change  in  T2  with  t  already  described.  That  is,  if 

the  spectra  of  deep  and  shallow  events  have  the  same  shape  and  see  the  same 

frequency  dependent  Q  along  the  raypath,  then  the  change  in  T2^S)  is 

* 

simply  that  required  to  account  for  the  change  in  t  .  Tentatively, 
this  result  is  interpreted  as  support  for  the  basic  similarity  between 
deep  and  shallow  shear  body-wave  spectra.  A  reasonable  corollary  is  that 
P-wave  spectra  are  also  similar. 

Under  the  assumption  of  similarity,  the  increase  in  i2(P)  with  source 
depth  must  be  interpreted  as  a  mixed  depth  and  frequency  dependence  in  Q^. 


Specifically,  an  isolated  bulk  loss  mechanism  must  operate  in  the  upper 


mantle.  The  fact  that  T^(P)  increases  toward  t2(S)  implies  that  the 
amount  of  bulk  loss  seen  along  the  travel  path  is  decreased.  An  intuitively 
appealing  explanation  is  that  the  bulk  loss  mechanism  operates  within  the 
asthenosphere  where  a  partial  melt  is  hypothesized. 

These  conclusions  may  be  visualized  in  Figure  3.5.  The  solid  line 
represents  the  right  side  of  the  shear  absorption  band  as  it  affects  P-waves. 

The  dotted  line  represents  the  observed  P  wave  absorption  band  which  is  a 
composite  of  bulk  and  shear  losses.  The  dashed  line  represents  the  hypothetical 
bulk  loss  mechanism  which  controls  TjCP)  for  shallow  focus  events.  For 
depth  of  focus  beneath  the  bulk  loss  region,  only  half  as  much  bulk  loss 
is  seen  by  a  ray,  effectively  pushing  T2(P)  toward  the  limit  of  the  shear 
absorption  band.  Note  that  the  implied  relaxation  time  for  the  bulk-loss 
mechanism  as  sketched  is  about  0.5  sec. 

Frequency  Dependent  Earth  Models. 

The  results  of  the  preceding  section  provide  a  starting  place  in  the 
construction  of  a  frequency  dependent  Q  .  In  the  examples  discussed 
below,  T2(P)  will  be  assumed  to  be  controlled  by  shear  mechanisms  in  the 
crust  and  in  the  mantle  below  the  asthenosphere.  In  the  asthenosphere, 
where  a  partial  melt  may  exist,  is  assumed  to  be  controlled  by  a 

bulk  loss  mechanism.  The  models  AFL  and  SL1  will  provide  a  framework  for 
model  evaluation.  As  an  example  of  the  flexibility  offered  by  manipulation 
of  T2  ,  preliminary  attempts  will  be  made  to  generate  AFL  from  SL1  by 
adjusting  the  position  and  strength  of  a  bulk  loss  mechanism  in  the 
asthenosphere. 

* 

The  difference  in  average  t  levels  between  AFL  and  SL1  is  easily 
modeled  by  adjusting  T2  without  depth  dependence.  Table  3.5  gives  four 


P-WAVE  ABSORPTION  BANDS 


Figure  3.5:  t*  vs  Frequency.  The  effect  of  an  isolated  bulk  loss  mechanism  is  shown.  For  shallow  focus  events, 
the  bulk  loss  makes  the  P-wave  absorption  band  wider  than  the  shear  loss  absorption  band.  P  waves 
from  deep  events  see  only  half  as  much  bulk  loss,  so  that  the  P-wave  and  S-wave  absorption  bands  are 
more  nearly  the  same. 
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TABLE  3.5 

FREQUENCY  DEPENDENT  Q 
MODEL  PARAMETERS  AND  DISTANCE  DEPENDENCE 


f 

! 

1 

Model 

1 

Model 

2 

Model 

3 

Model 

4 

Depth 

T2(P) 

Depth 

T2(P) 

Depth 

t2(p) 

Depth 

t2(P) 

1 

0 

0.26 

0 

0.18 

0 

0.26 

0 

0.18 

i 

2886 

0.26 

61 

0.18 

61 

0.26 

61 

0.12 

t 

j 

\ 

t 

81 

0.08 

81 

0.16 

81 

0.09 

196 

0.08 

196 

0.16 

196 

0.09 

f 

221 

0.18 

221 

0.26 

221 

0.10 

2886 

0.18 

2886 

0.26 

321 

0.12 

371 

0.14 

388 

0.18 

2886 

0.18 

A(DEG) 

* 

fccx 

200  sec 

1  sec 

* 

ca 

200  sec 

1  sec 

* 

200  sec 

1  sec 

* 

’a 

200  sec 

1  sec 

80 

.800 

.284 

.771 

.400 

.770 

.274 

.772 

.415 

70 

.858 

.302 

.822 

.426 

.822 

.297 

.824 

.442 

60 

.984 

.333 

.913 

.470 

.912 

.318 

.914 

.487 

50 

.942 

.327 

.918 

.472 

.917 

.350 

.919 

.490 

40 

.860 

.304 

.868 

.451 

.865 

.333 

.868 

.472 

30 

.755 

.267 

.761 

.405 

.760 

.302 

.762 

.431 

Note:  t  values  are  for  depth  of  focus  =  5  km. 
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examples  of  vs  depth  which  have  been  superimposed  on  SL1  and  the 

j k 

corresponding  t  values  for  waves  of  period  200  sec  and  1  sec.  From 

Figure  3.5,  it  is  obvious  that,  for  frequencies  greater  than  w  =  I/t^  , 

-2  -3 

increasing  T2  increases  the  Q  .  Both  12(0)  )  and  t2(q)  )  give 

* 

approximately  correct  change  in  t  with  frequency. 

* 

The  difference  in  t  vs  hypocentral  depth  is  not  adequately  modeled 
by  any  of  the  distributions  shown  in  Table  3.5.  t*  vs  depth  is  shown 

in  Table  3.6. 

The  complete  inversion  for  (depth)  required  to  obtain  model  AFL  from 

model  SL1  will  be  left  to  a  future  report.  But  it  is  apparent  from  the  testing 
already  done  that  the  result  will  be  almost  trivial.  At  each  depth 

Q_1(AFL)  =  Q_1(SL1)  •  R(u,  T2) 

Q(AFL)  =  1 

Q(SL1)  R(w,  T2) 

Note  that  the  ratio  of  high  frequency  Q  to  low  frequency  Q  must  be 
greater  than  or  equal  to  one,  unless  a  bulk  loss  is  operating.  For  these 
two  models,  the  ratio  is  less  than  one  for  a  depth  range  0-160  km,  and  is 
always  greater  than  two  below  160  km.  Again,  the  existence  of  a  bulk  loss 
concentrated  in  the  asthenosphere  is  supported  by  the  data. 

Summary 

A  review  of  published  Q  models  was  undertaken  in  an  attempt  to  reduce 
the  number  of  variables  in  a  study  of  the  effect  of  frequency  dependent 
Q  upon  body-wave  spectra.  However,  rather  than  finding  a  consensus  among 
Q  models,  a  frequency  dependence  was  found,  corresponding  to  the  data  set 

it 

used.  For  the  purposes  of  further  work,  t  will  be  assumed  to  follow  the 


I 
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TABLE  3.6 

FREQUENCY  DEPENDENT  Q 

SOURCE  DEPTH  DEPENDENCE 

(%  Surface  Focus  t*) 
a 


Model 

1 

Model  2 

Model 

3 

Model 

4 

Depth 

200  sec 

1  sec 

200  sec 

1  sec 

200  sec 

1  sec 

200  sec 

1  sec 

5 

100 

100 

100 

100 

100 

100 

100 

100 

50 

97.9 

100 

99.6 

99.6 

99.9 

100 

99.6 

99.8 

100 

94.6 

97.2 

96.1 

95.1 

96.1 

95.7 

95.6 

95.3 

150 

91.0 

93.2 

92.4 

89.8 

92.5 

90.5 

92.3 

90.4 

200 

88.2 

90.4 

89.8 

86.2 

89.7 

87.1 

89.7 

86.9 

300 

83.5 

85.7 

85.3 

82.2 

85.2 

83.1 

85.2 

81.4 

400 

80.7 

82.3 

81.8 

78.8 

81.8 

80.0 

81.7 

77.8 

600 

76.7 

78.2 

77.2 

74.8 

77.3 

75.7 

77.3 

73.9 

Note:  t  values  are  for  A  =  50° 
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f ree-oscillation  Q  model,  SL1.  This  choice  is  based  upon  the  low-frequency 
data  used  to  compute  SL1  and  upon  the  fact  that  the  frequency  dependence 
used  in  this  study  includes  no  variation  in  Q  at  low  frequencies. 

Given  a  starting  model,  then,  the  approximation  of  other  Q  models  by 
adjustment  of  vs  depth  is  easy.  Indeed,  any  high-frequency  Q  model 

may  be  generated  from  any  low-frequency  model  by  adjusting  the  width  of  the 
absorption  band.  Rather  than  a  single  infinity  of  Q  models  which  will  satisfy 
a  given  data  set,  there  is  now  a  double  infinity  of  models,  even  though  the 
frequency  dependence  modeled  is  simple  enough  to  be  characterized  by  a  single 
parameter . 

It  is  thus  even  more  important  to  constrain  Q(f,  r)  with  new  data. 

Of  particular  importance  are  the  body-wave  spectra,  where  direct  observation  of 
may  be  made.  Time-domain  modeling  of  pulse  shapes  may  also  improve  the 
constraint.  But  basically  any  ray  samples  an  average  of  properties  along 
the  path,  and  detailed  inversion  for  depth  dependent  properties  cannot  be 
done  simply.  The  conclusion  put  forth  in  this  paper  is  that  only  the  end  of 
a  complicated  band  of  relaxation  mechanisms  may  be  resolved  by  body-wave 
spectra.  If  the  window  in  the  spectrum  of  relaxation  mechanisms  were  not  in 

the  passband  of  stands  \  seismometers,  then  not  even  the  end  of  the  absorption 
band  could  be  observed. 
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IV.  Applications  of  Green's  Function  Methods  in  Elastodynamics  to  Source 
Theory  and  Wave  Propagation  Problems 

C.B.  Archambeau 

A  major  part  of  the  effort  during  the  past  six  months  of  the  contract 
period  has  been  to  develop  a  comprehensive  theoretical  framework  for  the 
treatment  of  both  source  and  wave  propagation  problems  in  elastodynamics. 

The  Appendix  1  provides  a  summary  of  the  Green's  fucntion  theory  which  will 
provide  such  a  framework  for  elastodynamic  problems  in  Geophysics.  In 
this  theoretical  development  we  consider  media  which  nay  be  inhomogeneous, 
anisotropic  and  have  both  fixed  and  moving  discontinuities  in  material 
properties,  the  latter  being  included  in  order  that  generalized  relaxation 
sources  in  prestressed  media  may  be  treated.  This  work  corresponds  to  a 
generalization  of  previous  work  (e.g.,  Archambeau  and  Minster,  1977) 
mentioned  in  earlier  reports,  which  was  specifically  focussed  on  source 
problems.  The  present  work  is  more  general  in  that  we  have  included  both 
fixed  and  moving  boundaries  in  the  same  formulation  and  most  importantly, 
have  generated  explicit  Green's  functions  for  layered  media  to  be  used  in 
the  integral  equations  for  applications.  We  are  therefore  in  a  position  to 
solve  a  large  number  of  outstanding  problems  in  seismic  source  and  wave 
propagation  theory  and  specifically  particular  problems  of  importance  for 
seismic  event  discrimination  and  explosion  yeild  estimations. 

In  this  regard,  in  this  section  we  show  first  how  layered  half  space 
Green's  functions  are  obtained  and  give  explicit  relations  for  them. 

We  then  show  how  such  a  Green's  function  (or  more  properly  this  layered 
half  space  Green’s  tensor)  can  be  used  with  the  integral  formulations  in  the 
Appendix  2  to  provide  the  means  of  representing  a  general  nonlinear  (or  linear) 
energy  source  imbedded  in  a  layered  half  space  (which  may  be  specified 


numerically  in  general)  in  analytic  form.  Such  a  result  can  then  be  used 
to  predict  both  surface  and  body  wave  radiation  from  general  seismic  energy 
sources;  in  particular  from  complex  numerical  models  of  both  earthquakes 
and  explosions. 

In  addition,  we  have  formulated  the  same  theory  in  a  layered  spherical 
medium.  Since  this  theory  is  very  similar  in  principal  to  the  layered 
half  space  theory  we  do  not  include  the  details  here. 

Other  applications  of  this  theoretical  formulation  are  being  studied, 
in  particular  effects  of  lateral  variations  in  iructure  on  predicted 
source  radiation,  effects  of  scattering  from  the  growing  source  boundary 
and  a  variety  of  perturbation  modifications  to  be  used  to  predict  mroe 
accurately  the  effects  of  anelasticity .  velocity  gradients  and  boundary 
variations  from  spherical  or  planar  form.  These  additional  applications 
will  be  described  in  subsequent  reports. 
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Eigenvectors  for  Layered,  Isotropic  Elastic  Half-Spaces 

Eigenvectors  corresponding  to  the  operator  defined  in  the  Appendix  II, 
equations  (61)  and  (62)  are  generated  by 


Vk +  "A 


(4.1) 


In  a  layered  space,  with  boundaries  along  the  coordinate  surfaces 
**  zs»  s  »  0,...,n-l,  the  set+ 

n 

1 


L(3)  + 


(s) 

k 


(s)  2  (s) 

p  to  ip 


) 


(4.2) 


is  equivalent  to  (63).  That  is  is  represented  by  the  differential 
operator  set  j  j  and  iji^  by  the  eigenvector  set  {^S^}  *  The  members 


4-Indices  enclosed  by  parentheses  are  not  subject  to  the  summation  rule.  Thus 
no  sum  is  implied  on  the  repeated  indices  (s)  for  example. 


28 


of  the  sets  are  connected  by  the  boundary  conditions,  which  are,  from  (62): 


z 

o 


Here  the  free  surface  is  the  coordinate  plane  x_  =  z  .  These  conditions 

3  o 

insure  that  the  members  of  the  eigenvector  set  produce  a  continuous  eigenvector 
^  with  associated  continuous  tractions.  Explicitly,  ^  is  given  by 


(s) 

with  the  adjoining  continuous  at  the  boundary  points  z  .  A  similar  statement 

applies  to  the  tractions. 

Now  the  elastic  parameters  and  density  are  constant  within  each  layer 
by  definition,  so  that  requiring  isotropy  as  well  produces  the  result: 


a2  *(s)  . 


(4.4) 


;  s  =  0,1,. . ,n-l 
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where 


B 


(s) 

ilk 


(1_<Sso)5jk- 


with  8^  defined  by 


(A.  5) 


B(s),p(s) 
Jlk  vk 


+  y 


<s)  f  Jjj _ +  — 

**t  Sxj 


(A. 6) 


which  is  just  the  traction  '^|C£n£*  Here  the  external  boundary  condition  has 
been  incorporated  with  the  boundary  conditions  on  the  internal  boundaries 
by  defining  A^  =  y^  =  p ^  =  0.  In  addition  the  6x3  matrix  operator 
has  been  introduced  for  convenience  in  writing  the  boundary  conditions 
in  compact  form. 

The  vector  form  of 


+  y*s^]  v(v-i(s))  +pS^  i(s)  +p(s¥s)=0 


(A. 6a) 
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From  (A. 4)  or  (4.4a)  it  is  easy  to  see  that  solutions  for  may  be 
generated  by  introduction  of  the  (physical)  potentials  defined  by 


1 

2  Ej  kl 


(4.7) 


So  that  these  potentials  can  be  represented  by  a  four  vector  with  cartesian 
(s) 

components  X^aj>  a  a  1,2, 3, 4.  In  vector  form: 


1/2  Vx 


(4.7a) 


Using  these  in  (4.4)  or  (4.4a)  shows  that  if  the  cartesian  components 

(s) 

satisfy  the  scaler  Helmholtz  equation: 


y2Y(s)  + 
X(«) 


a  -  1,2, 3,4 
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then  ^  is  given  by 


■  (kis).  kis).  *.<a).  -<a>) 


where  k' 


’  &mn  Sx 


(4.8) 


o»/v(s)  ;  v(s) 

P  P 


+  2u(s* 


*»  u/v^ 
s  s 


00  .  00 
„  »  V 

B  8 


/P(S)/P(S) 


Alternately,  In  vector  form,  /s)  is  generated  from  X<S>  and 
X  ■  (X1,X2,X3)  by 


■[kP<8>]  %  +  2[k,<s)] 


(4.8a) 


provided  the  potentials  x4  and  x  satisfy  the  scalar  and  vector  Helmholtz 
equations  respectively.  That  is,  when 


v2  00  + 
x4  + 


(4.9) 
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I 


Eigenfunctions  for  the  cartesian  components  of  the  potentials  are 
the  set 


„<s> „ 

(kp)  elm* I  e  M 


J  (kp)  elB* 
m 


and 


jvkp)  *lm*j k+  2 


where 


(s)  =  ^k2  -  [k<s)J2 


00 


» 


It  is  critical  to  this  development  that  it  turns  out  that  only  depends 

on  the  layer  parameters,  as  evidenced  by  the  index(s). 


(s) 

Therefore  the  cartesian  components  of  can  be  written  as: 


00 

(°0 


v<8)z 

A(S)  (k  .)  e  <a) 


_v(s)2 

+  bJs>  (k,u)  e 
(a).m' 


J  (kp)  e 
m 


im<p 


(4. 


10) 
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The  reason  for  the  introduction  of  coefficients  in  the  expressions 

for  the  Y  is,  of  course,  that  iji  must  satisfy  boundary  conditions  in  order 
a 

that  it  be  an  eigenvector  of  the  operator  defined  by  (63)  and  (64),  and 
the  coefficients  are  to  be  used  to  satisfy  these  conditions. 

f  g  \ 

Use  of  (4.10)  in  (4.8)  produces  the  cartesian  components  of  ijr  .  These 

can  be  expressed  as  components  in  a  cylindrical  basis,  which  is  most 

convenient  here.  In  particular,  the  eigenvectors  can  be  expressed  in  terms 

of  the  vector  cylindrical  harmonics  P  ,  B  and  C  as  (Ben  Menahem  and 

— m  - m  — m 

Singh,  1972j  Morse  and  Feshbach,  1953); 


^s)(r,k)  -  D‘S)(Z.W  T,+  4SW>  *„  +  C, 


(4.11) 


where 


P  -  €  J  (kp)e 
— m  z  m 


im<j) 


-  e 


■  ^ 

- +  e  . 

P  3(kp)  * 


c 

—TO 


t 

[ 6p  (k) 


(i)  I*]  V 

K 


9_  „  3 _ 

3$  "  e4>  3(kp) 


im$ 


im<{> 


(4.12) 


Here  P  «B  -  P  *C  -  B  »C  -  0. 
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The  coefficient  functions  of  z  are  (Ben-Menahem  and  Singh,  1972) 


(s)  (s) 

D<s)(z,k)  =  v(s)  -a(s)(k)  e  p  +  b(s)(k)  e  p 
m  p  m  m 


„(s)  ..(s) 

+  k  c<s><k>  e  s  +d(s)(k)es 
in  m 


, ,  , ,  -JSK  , ,  JSK 

E(s)<z,k)  -  k  a<s)(k)  e  P  +b<s)(k)e  p 
ram  m 


v(s)  (s) 

-c(s)<k)e  8  +  d(s)(k)  e  8 

m  m 


,  .  ,  .  -v<s)z  ,  .  v(s)z 

F(s)(z,k)  -  «iS)Ck)  e  +  !() (k)  e 

m  m  ra 


(s)  (s) 

The  coefficients  a  ,  b  ,  etc.  are  linear  combinations  of  the  coefficients 
m  m 


(4.13) 


appearing  in  the  expressions  for  the  potentials 
Note  that: 


k2  -(k<*>)2  ;  k  >  k<*> 

^  ‘  ^  '  kp 


and  similarly  for  .  Also,  the  radiation  condition  requires  that: 
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b(n)  =  d(n)  =  f(n)  =  0. 
m  m  m 


The  boundary  conditions  of  (4.3)  must  be  satisfied  in  order  that  (4.11) 
be  an  eigenvector.  Using  (4.11)  in  these  conditions  shows: 

(1)  that  the  vector  function  P  ,  B  ,  C  can  be  eliminted  from  the 

— m  -m  — m 

boundary  condition  equations,  so  that  only  the  coefficients 
(s) 

D  (k,z) ,  etc.  are  involved  in  the  boundary  conditions, 
m 

(2)  The  boundary  equations  at  each  layer  interface  can  be  separated 

into  two  independent  sets,  the  first  conyosed  of  the  coefficients 

of  P  and  B  in  the  second  involving  only  the  coefficient 

— m  — m 

of  Cm»  These  correspond  to  the  P-SV  coupled  waves  and  the  SH 
waves  respectively.  They  therefore  also  represent  Rayleigh  and 
Love  type  surface  wave  modes  when  the  eigenvectors  are  viewed  as 
inodes  of  oscillation  of  the  medium. 

as  the  sum  of  two  decoupled  eigenvectors 
Thus 


The  eigenvector  can  therefore  be  viewed 
which  may  be  denoted  by  *S|/S^and  . 


R^s)  +  Ljj/S> 


=  D^P  +  E^B 
m  ~m  m  ~m 


VS)=F(S) 

•tb  m 


C 

— m 


(4.14) 
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All  the  boundary  conditions  are  now  expressed  by  the  P-SV  and  SH  boundary 
conditions.  In  matrix  form  they  are: 

(1)  P-SV  boundary  conditions: 


D(s)(z) 
m 


E(s)(z) 

m 


p(s) (z) 
in 


E(s)(z) 

in 


d(8+d  (z) 

m 


E(s+1)(z) 

m 


P(s+1)(z) 

m 


E(s+1> (z) 
m 


(4.17) 


(2)  SH  boundary  conditions 


FiS) (z) 

m 


FmS)  (2) 
m 


F(s+1>(z) 

m 


F(s+1)  (z) 
m 


(4.18) 


The  equations  involving  the  unknown  coefficients  a^  ,  b^  »•••  etc.  in  (4.13) 
and  (4.16)  may  also  be  written  in  matrix  form.  In  particular,  at  z  =  zg 
we  have,  for  the  P-SV  coefficients: 


DiS)(z) 

m 


Eis)  <z> 

m 


» 


m 


KpS>(*  ) 
R  s 


(4.19) 


,(s) 


m 


I 


J 
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and,  for  the  SH  terms,  similarly: 


z 

s 


(s)  (s) 

The  matrices  and  are  4x4  and  2x2  respectively,  and  are 

invertable  (e.g.  Haskell,  1953;  Harkrider,  1964).  Here  of  course  the 
(s  )  (s  ) 

coefficents  a  ...f  are  not  functions  of  z.  Therefore  we  can  evaluate 
m  m 

the  equations  at  z  =  z  ,  where  this  is  the  upper  boundary  of  the  layer 
(s)  and  get  another  expression  for  the  constant  coefficients.  Thus 
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(s)  (s) 

These  may  be  solved  for  the  coefficients  a  ...f  in  terms  of  the 

m  m 

functions  of  z  ,  and  one  has: 
s-1 


D^(z 

m  '  s-1 


E^s) (z  ,) 

m  s— 1 


(4.23) 


(4.24) 


Now  we  can  eliminate  the  coefficients  from  (4.19)  and  (4.20),  this  giving  us  a 
result  which  effectively  propagates  the  solution  across  the  layers,  from 
the  zg  j  boundary  to  the  zg  boundary.  In  particular,  from 


‘Vl* 


(4.25) 
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and  from  (4.20)  using  (4.24): 


where  these  are  the  Rayleigh  and  Love  or  P-SV  and  SH  propagators. 


where  the  notation  ha9  been  modified  in  an  obvious  way  for  brevity 


(4.26) 


(4.27) 


Then 


(4.28) 


Now  the  boundary  conditions  require 


the  final  equality  because  of  (4.28).  ^Using  (4.28)  again,  with  (s)  taken  as 
(s-1)  in  the  expression, gives 


for  example,  with  a  similar  expression  for  the  SH  propagator  equation. 
This  can  be  used  in  (4.29)  to  eliminate  the  factor 


and  likewise  for  the  SH  equation.  Clearly  this  can  be  continued  until 
we  reach  the  z  boundary  at  the  free  surface.  We  get  in  fact: 
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(4.30) 


The  final  boundary  condition  to  be  met  is  the  vanishing  of  tractions  at  the 

free  surface  z  =  z  =  0,  so  0^(0)  =  E^  (0)  =  (0)  =  0  in  (4.4.30).  On 

o  m  m  id 

the  other  hand  the  displacements  are  non-zero  at  z  -  0.  One  of  them  in  the 
P-SV  eigenvector  may  be  set  to  unity  for  the  eigenvector  solution  sought 
here,  with  the  eigenvector  normalization  factors  used  to  account  for  this 
later.  For  the  decoupled  SH  vector,  there  is  only  one  displacement 
amplitude  and  so  it  may  also  be  set  to  one, with  normalization  accounting 
for  it  later.  Thus  we  have 
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with  e  =  (0)  (0) .  (Use  of  e  as  defined  amounts  t  o  dividing  both 

o  m  m  o 

sides  of  (4.30)  by  (0)  for  P-SV  and  by  F^(0)  for  SH,  and  then 

m  in 

redefining  the  coefficients  so  that  D’  a  ^  (0)  •  Because  of 

m  mm 

our  freedom  to  normalize  the  eigenvectors  however, we  can  simply  write 
these  new  coefficients  as  before.) 

Therefore,  if  this  relation  is  satisfied  for  s  =  n-1,  then  all  the 
boundary  conditions  are  satisfied.  (Mote  that  [R^  ^..R^  ^ ]  must  be  taken 
as  unity  when  s  =  0.) 

The  coefficients  required  to  insure  that  (4.31)  is  satisfied  for  all 
s  «  0,...,n-l  are  obtained  from  (4.23-4.24)  using  (4.31).  That  is,  we  require 
that  (4.32)  be  valid  for  any  s,  in  particular  for  s  =  p-1.  The  same 
holds  for  (4.23)  and  (4.24)  for  s  «  p.  Therefore  taking  s  =  p  -  1  in  (4.31) 
and  substituting  the  resulting  relation  in  (4.23)  and  (4.24)  with  s  =  p 
used  in  these,  we  have 


m 


when  p  >  2. 
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Defining: 


Here  Jp  is  a  4  x  4  matrix  propagator  and  JjP^  a  2  x  2  matrix  propagator. 

K  •4 


(*lo 

p  o 


J(p)  H  J(p)  (z  Iz  > 
L  L  P  ° 


That  is: 
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propagate  the  solution  vector  from  the  boundary  surface  at  2  to  the 

o 

boundary  at  z^. 

Explicitly,  the  coefficients  are: 


(4.35) 


(■'■').(  [-’’'I) 

v-y  Ip:"],./ 

With  the  coefficients  satisfying  these  relations,  then  the  solutions  in 
(4.13)  are  in  fact  eigenvectors  of  the  operator  of  (4.1)  and  (4.2). 
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Here  the  original  differential  equation  has  been  divided  through  by  p, 
the  density,  in  order  to  provide  the  same  differential  operator, 
p  that  was  actually  used  to  generate  the  eigenvectors  of  (4.37). 

The  functional  orthogonality  of  the  i(j(i:,k)  can  now  be  easily  demonstrated 
and  this  in  turn  can  be  used  to  obtain  the  expansion  for  H™  in  terms  of 
these  eigenvectors.  (Actually,  completeness  of  the  eigenvector  set  is 
also  needed  for  such  an  expansion  to  be  rigorous.  However,  we  will  obtain 
the  expansion  without  the  proof  of  completeness,  just  using  the  orthogonality 
of  the  Since  the  it™  so  generated  satisfy  both  the  differential  equation 
for  a  Greens  function  and  the  required  boundary  conditions  then  by 
uniqueness  of  solutions  for  such  an  equation  we  can  claim  to  have  a 
sufficiently  proper  representation  in  the  distributional  sense.) 

The  orthogonality  of  the  eigenvectors  viewed  as  a  set  composed  of 
members  with  different  eigenvalues  m  and  k,  is  a  consequence  of  the  self¬ 
adjointness  of  the  operator  defined  by  (all)  the  equations  of  (4.39). 

(Here  m  and  k  are  eigenvalues  or  "quantum  numbers"  for  eigenvectors  of 
the  operator.  Clearly  m  is  a  discrete  set  of  integers,  while  k 
is  a  continuous  set.  For  a  bounded  medium,  like  a  layered  sphere,  the 
eigenvalues  appearing  are  m  and  l.  They  correspond  (one  to  one)  to  the 
set  m,  k  appearing  in  this  treatment  of  an  unbounded  medium.  For  the 
bounded  spherical  medium  however  both  m  and  1  are  discrete  integer  sets. 

These  differences  are,  as  implied,  a  consequence  of  the  bounded  or  unbounded 
nature  of  the  space  over  which  the  operator  is  defined,  the  basic  operator 
being  the  same  in  both  cases.) 

The  self-adjointness  of  (4.39)  is  demonstrated  by  considering  two 
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arbitrary,  distinct  members  of  the  eigenvector  set  j}>  and  (say)  _£>  where 
we  use  different  symbols  for  these  members  for  clarity.  Since  and  £ 
are  distinct  they  will  have  two  distinct  sets  of  eigenvalues.  In  particular, 
the  eigenvectors  satisfy: 


Lik*k 


(4.40) 


where  L’  *»  p  is  used,  with  the  distinction  of  the  prime  ignored. 

2  2 

Here  both  w  and  v  are  real.  The  inner  product  of  with  the  first  of 
these  gives  (see  Archambeau  and  Minster,  1977  for  details): 


where 


(4.41) 


Jk  “  **  [Ciklj8J*i  ]  “  [C»kij8j*i] 


(4.42) 


is  the  bilinear  form.  Here  denotes  the  complex  conjugate  of 

is  the  adjoint  differential  operator  and  is  found,  in  Appendix  II,  (49)  to  be  a 
self-adjoint  differential  expression,  in  that: 


(4.43) 
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The  inner  products  are  defined,  for  example,  by 

(Wfh)  -f 

Now  integrating  the  last  term  using  the  divergence  theorem,  taking 
account  of  the  existence  of  internal  boundaries  in  v,  gives: 


+f  L £  (8£i'l'i)  “  da 

3vz  ■* 

where  is  the  boundary  operator  defined  in  equation  (68) .  Now  and  <p 
are  eigenvectors  satisfying  the  boundary  conditions  of  (102), and  (100) 
plus  (101)  is  their  explicit  form.  Thus  since 
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then  the  integral  over  the  internal  boundaries  vanishes.  Likewise,  since 


then  the  first  integral  vanishes  as  well.  We  therefore  have  that 


/Jk,kd3*  - 0 


We  have  as  a  consequence  that: 


Since  ^  and  satisfy  (103),  then  this  gives 

2  2  *f 

But  w  4  v  ,  since  and  £  are  distinct,  so  that  we  must  have  : 


40n  occasion  we  will  use  to  denote  inner  products  as  well  as  (4*,+}. 


(4.44) 


•/ 

v  J  v 


Vi 


d3x 


0 


Hence  and  £  are  orthogonal  in  both  the  vector  inner  product  sense  and 

the  functional  sense,  (i.e.,  the  repeated  index  &  requires  summation 

over  its  range  so  that  we  also  have  a  vector  inner  product  on  v. ) 

The  vectors  of  the  set  are  therefore  mutually  orthogonal,  however 

we  do  not  assume  that  they  are  normalized.  The  normalization  factors 
P  8  C 

N  (k,u),  N  (k,w),  and  N  (k,w)  may  be  functions  of  the  eigenvalues  m,  k 
mm  m 

and  are  associated  with  the  vectors  P  ,  B  and  C  of 

— m  — m  — m 

defined  from  the  magnitude  of  the  inner  product: 

<  s/’w3* 

In  view  of  the  vector  orthogonality  of  P  .  B  and  C  ,  we  can  define  norms 

—m  — m  — ra 

for  each  of  the  component  vectors  making  up  in  particular: 


HP(k)  -  /  P  •?  d3x 
m  f  - m  -to 

v 

NB(k)  =  f 8  •B  d3x 
m  I  ~ m  -m 

N^(k)  =  f  Cm’Snd3x 

in  f  in  -m 


\ 


(4.45) 
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In  view  of  the  fact  that  ip  is  the  sum  of  P  ,8  and  C  defined  in  (4.37), 

—  m  m  m 


then  the  orthogonality  of  £  is  expressed  by  (using  (4.44)  and  (4.45)  together) 


/?  -P*  ,d3x  =  NP6m,6(k-k’ 

— m  — m  mm 

v 

/  8  ‘B'.d^  =  H8Sm’6(k-k 

J  — m  — m  m  m 

/C  *C\d3x  =  NC6m'6(k-k') 
— m  — m  m  m 


*) 


(4.46) 


The  orthogonality  of  the  eigenvectors  can  be  used  to  provide  an  explicit 
expression  for  the  Greens  tensor  Hq.  In  particular  we  have  that 


satisfies  an  equation  of  the  type 


L  Hq  +  v2^  = 
rp  p  r 


-  4ir646(r-r  ) 
r - o 


Here  v  is  the  Fourier  time  transform  parameter,  corresponding  to  the  angular 


frequency.  Also  satisfies  the  boundary  condition  of  the  form 


We  take  H^(r,r  )  to  be  expressible  by  an  expansion  in  the  eigenvector 
p - o 


set  <|Kr) .  For  brevity  we  represent  the  eigenvalue  pair  m  and  k  by  the 
single  symbol  X.  Further  we  write  the  sum  over  the  set  of  eigenvectors  as 
a  "sum"  over  the  sets  of  m  and  k,  represented  by  /  .  Since  m  is  integer 


and  discrete  while  k  is  a  continuous  set  then  /  is  equivalent,  in  this  case, 

X 


to  an  ordinary  sum  over  the  discrete  set  of  m  values  after  integration  over 
the  continuous  k  set.  Thus  we  take: 


■ 
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^(r;r  ) 
p - o 


Aq(VX)  *p^jX) 


A 


(A. 47) 


where  the  A(r  ;A)  are  unknown  coefficients  for  the  expansion.  They  may, 

— o 

of  course,  be  functions  of  r  and  A  as  indicated. 

*  — o 

Substituting  this  expansion  in  the  differential  expression  for  gives 


J  »x)  I" +  =  ~  5(r-r  ) 

./  q  o  rprp  rl  r  - o 

\ 

Taking  the  inner  product  with  another  member  of  the  eigenvector  set 
and  observing  that 


L 


rp  p 


2. 
id  ip 


r 


in  the  above,  we  have 


/VVA)  j^((02-v2)  <  ij»r  |  <f>r  >  j  =  An  <«5U-ro)  | 


<frq(r;A')  ^ 


(4.48) 


Using  this  and  the  orthogonality  of  the  eigenvectors,  we  can  determine  the 
coefficients  of  A  (r  ;A)  in  the  Green's  function  expansion  (4.47). 

q  — o 

However,  to  obtain  normalization  factors  for  the  A  in  the  Green's 

q 

function  expansion  which  are  consistent  with  the  use  of  the  ratio 

(ellipticity  factor  )  eQ  in  reducing  the  eigenvector  coefficients  in 

equation  (4.30)  to  those  of  (4.31),  we  must  express  the  eigenvectors  as  the 

R  L 

sum  of  a  P-SV  type  vector  and  a  SH  vector,  in  the  form  ,  and 
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R  L 

normalize  ^  and  ^  separately.  This  decomposition  simply  means  that  the 

2 

equations  of  motion  £.  <»  +  w  <J»r  =  0  and  boundary  conditions  can 

simultaneously  be  broken  into  two  independent  sets  with  one  set  satisfied 
R  L 

by  £  ,  the  other  by  £  .  The  analysis  starting  with  equation  (4.37)  and 
resulting  in  equation(4. 38)  applies  to  each  operator  set  (differential 
equation  plus  boundary  conditions)  involving  and  .  Then  H™  is  given 
by 


u®  =  Hm  + 

Hi  RH  i+  Li 


with 


R2£k 


R^k  + 


2  rim 

ui  =  - 

R  i 


4pd™6R(r  - 


r  ) 


2  rim 

*  iBi 


-  -  I„> 


where  6(r-r  )  =  6„(r-r  )  +  6_ (r-r  ),  in  which  6  and  6  are  decompositions 
of  6  such  that  <dR!j^L>  =  0  and  < <5 ^ | =  0. 

We  are  then  led  to  results  parallel  in  all  respects  to  those  already 
obtained  for  the  eivenvectors ;  and  such  that 
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with  6R(_r  ~  £q)  and  6^0:  -  £q)  being  regular  delta  functions  of  jr  and 

r^,  but  with  the  added  property  that  these  are  projection  operators  along 
R  L 

and  £  (e.g.,  Morse  and  Feshbach,  1953).  Here: 
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Further 


<«R(r-ro) 


<frq(r;X')  >  =>  ^(r^X') 


<6.(r-r)  *L(r;X')  >  =  ^(r  ;X’) 

L>  —  — o  q  —  q  — o 


Therefore  ( 111-a. )  and  (  111-b)  give 


AL(r  ;X')  =  —  -  —  •  ;X’) 

q  -O  NR(X«)  * 


vvx’> 


4L(r  ;X’) 


nV)  q~° 


Here  we  have  set: 


NR(X')  =  (/,+  N®,)C“2  ~  v2) 


m  m' 


NL(X’)  E  (/,x<*2  -v2) 


The  notation  $  indicates  the  complex  conjugate  of  <t>. 


(4.51) 
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R  L 

Since  >=  £  +  is  just  one  of  the  eigenvectors  of  the  set  and  X'  its 
associated  eigenvalue  set,  then  it  is  equivalent  to  write 


(4.52) 


The  Green's  function  is  therefore  given  by 


AR(r  ;X)  =  —  <t£(r  ;X) 

q  ~°  nR(x)  q  “° 


AL(r  ;X)  =  —  ^(r  ;X) 

q  “°  NL(X)  q  “° 


Hq(r;r  ) 

p  o 


_Hq  +  THq 
R  p  L  p 


"/ 


NR(X) 


^<VX> 

NL(X) 


(4.53) 


writing  this  out  more  fully,  using  the  appropriate  representation  of  the 
generalized  sum  over  X,  we  have 
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OO 

H’dsr,)  -  4,  £  /dkfkr[(V2=k> 

m  o  >  in'  L ' 


+  E  (z;k)  B 
m  *  —i 


D„(z„;k)  P  (p,i  ;k) *e 
o  tn  o  o  q 


(p^;k).;p^Dm( 

ik)\j] 

^<P,*:k>'lp)(?«<2o;k)‘  5.<Wk>'£<,)]} 


+  Em^2n»k^  it/Pjf’ 

in  o  m  o  o 


0 


I 

Here  ep  and  e^  are  the  p  and  q  components  of  the  coordinate  basis  chosen. 
Of  course  a  cylindrical  coordinate  basis  is  the  natural  choice  for  this 
representation. 

The  functions  D  ,  E  and  F  are  defined  in  equation  (4.11)  and  P  .  B 
mm  m  -m*  — m 

and  C  in  (4.10).  The  coefficients  in  the  function  D  ,  etc.,  are  given  by 
— m  m 

(4.35)  and  (4.36). 


.54) 


59 


Representation  of  the  Radiation  Field  from  a  Linear  or  Non-linear 
Energy  Source 

We  have  at  our  disposal  the  integral  representation  of  the  spectrum  of 
a  displacement  radiation  field  in  the  form  (equation  (60)): 

(A. 55) 


(A. 56) 


Consider  now  the  application  of  these  relations  to  the  representation 
of  the  elastic  radiation  from  a  general  energy  source  (e.g.  non-linear). 

We  will  assume  that  the  source  is  actually  given  by  a  detailed  numerical 
calculation  of  the  motion  due  to  some  particular  phenomena  (e.g.  an  explosion) 
which  is  most  likely  a  very  non-linear  process.  Suppose  the  source  region 
occupies  some  finite  volume  within  the  half  space  and,  for  generality,  that 
this  source  volume  or  nonlinear  zone,  may  intersect  the  free  surface.  The 
situation  may  be  as  shown  in  Figure  2.  In  any  case  v  is  the  linear  elastic 
zone,  which  will  be  assumed  to  be  a  layered  half  space.  The  external 
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and  the  layered  half  space  Green’s  tensor 


°  Air  y  .  I  dk 


00 

E/d 


tp(r;m,k)  ^(r^m.k) 

Ak) 

m 


*  (r;m,k)  if>  (r^jrn.k) 

^  T 

N“(k) 

tn 


60 


boundary  of  v,  denoted  as  3v  ,  can  be  viewed  as  being  composed  of  two 

L 

sections,  the  non-intersected  free  surface,  3vf  ^ ,  and  the  (elastic) 
boundary  between  the  (non-linear)  source  region  and  the  elastic  zone. 


With  the  layer  half-space  Green's  function  given  in  (4.56),  and  neglecting 
body  forces  at  least  for  the  moment (i.e.  the  first  integral  contribution 
can  always  be  added  later),  we  have  for  the  elastic  displacement  field  in 
v  (outside  the  source  region): 


4ttu  =  /  f  -  u  Hq  n  da 

q  J  L  p  pr  p  prJ  r  ° 


(4.57) 


3v 


(1) 


We  can  presume  knowledge  of  both  u  and  t  *n  ,  the  tractions,  on  3v 


(1) 


p  pr  r’  ’  E  * 

the  source  boundary,  since  in  the  present  application  we  have  supposed 

that  the  results  of  a  numerical  calculation  are  specified  on  3vj^  in  the 

form  of  Fourier  transforms  of  t  *n  and  u 

pr  r  p 

The  question  is  then,  how  does  this  excite  the  surrounding  elastic 
medium.  The  answer  of  course  is  given  by  (4.57)  in  the  form  of  the  surface 
integral  which  gives  u^,  the  elastic  displacement  at  any  point  in  v.  Ue 
note  that  the  integral  involves  only,  and  that  this  is  only  part 

of  the  nonlinear  zone  surface  boundary  that  does  not  intersect  the  free 
surface. 

The  representation  (4.57)  can  be  put  in  more  explicit  form.  That  is 
using  (4.56)  in  (4.57)  and  interchanging  the  integral  over  the  surface 


3v 


(1) 


with  the  sum  over  m  and  integral  over  k,  gives: 
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°  m  3v 
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* . 
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(4.58) 


*q(r;m,k)  f 

Ni<k)  i 

3v 


(u_)  C _ „  ~ ~ —  ^(r  ;m,k) 
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P  P^t  axt  Ts  -o’ 


n  da  > 
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We  observe  that  u  can,  as  expected,  be  written  as: 


u  -  uR  +  UL  (4.59) 

q  q  q 

Here  also  It  is  clear  that  we  can  define  tractions  associated  with  ^ 
the  eigenvectors  on  the  surface  Sv^.  (Previously,  in  equations  (4.13)  and 
(4.14),  tractions  on  planes  parallel  to  the  layer  interfaces  were 
defined  and  called  These  are  the  tractions  computed  for  the  layer 

matrix  and  used  in  the  propagators.  They  are  also  therefore  those  given 
by  the  standard  numerical  calculations,  e.g.  Harkrider,  1964.)  These  are 
undefined  until  we  choose  a  surface.  In  this  regard  we  have  some  flexibility 
of  choice  since  any  boundary  shape  so  long  as  it's  entirely  within  the 
linear  region  around  the  source  zone  will  do.  Clearly  a  surface  with 
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cylindrical  coordinate  planes  is  advantageous.  To  be  definite  about  the 

matter  we  will  assume  that  is  taken  to  be  a  cylindrical  surface. 

Thus,  the  eigenvector  stresses  when  denoted  as  ¥  and  ¥  ,  so  that: 

pr  pr’ 


H.R  «  .a  .  R 

4*  =  C  -  it 

pr  prst  axt  Ys 
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(4.60) 


YL  »  C  ~ 

pr  prst  dxt  s 
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can  be  used  to  provide  the  tractions  on  avf1^,  as  *n  and  4*^*  *n  .  These 

t.  pr  r  pr  r 

show  the  explicit  dependence  on  the  normal,  as  well  as  shorten 


(4.58).  That  is,  we  have 
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(4.62) 
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where 


-R  -L 

Therefore  u  and  u  have  the  forms: 
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where: 
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(4.63) 


These  coefficients  are  the  objects  of  interest,  since  once  they  have  been 

computed  it  is  a  rather  routine  matter  to  evaluate  (4.61a)  and  (4.62a) 

by  standard  procedures  —  such  as  by  summing  residues  in  the  k  plane  to 

obtain  surface  waves  —  and  by  evaluating  the  branch  line  integrals,  also 

occurring  in  the  k  plane,  by  approximate  methods  to  obtain  body  waves. 

(See  Ben-Menahem  and  Singh,  1972  for  examples.) 

(R  L) 

Since  we  know  the  forms  of  ^  *  ,  as  given  in  equations  (4,.ll)-(4. 13)  plus 

Cr  l) 

(4.35)  and  (4.36)  we  can  calculate  the  Ypr’  ‘n^  components  on  the  surface 
in  question  and  then  evaluate  the  integrals  in  (4.63)  by  a  suitable  numerical 
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integration  over  the  surface,  using  numerically  specified  fields  (u^)  and 

(t  *n  )  in  the  calculation, 
pr  r 


65 


Appendix  I 

Use  of  the  High  Order  Langer  Approximation  in  Seismogram  Synthesis 

V.F.  Cormier 

The  objective  of  this  research  is  to  develop  a  method  of  seismogram 
synthesis  that  is  more  efficient  in  computation  than  existing  methods  yet 
capable  of  predicting  waveform  modulation  by  zones  of  intense  vertical 
gradient  of  velocity  and  density.  In  routine  use  such  a  method  can  facilitate 
inversion  for  both  source-time  functions  and  earth  structure. 

The  last  semi-annual  report  considered  an  earth  having  regions  of 
intense  but  continuous  variations  in  elastic  properties.  Representations 
for  teleseismic  displacement  in  this  earth  were  derived  using  higher  order 
solutions  in  frequency  to  the  radial  eigenfunctions  that  solve  the  potential 
wave  equations.  Research  has  been  directed  now  towards  the  practical 
evaluation  of  displacement  using  this  representation.  The  higher-order 
solution  to  the  elastic  wave  equations  given  recently  by  Woodhouse  (1977) 
in  terms  of  propagator  matrices  is  shown  to  be  equivalent  to  the  solution 
derived  for  the  higher  order  potential  equations.  It  is  shown  that  for 
practical  calculations  it  is  best  to  describe  the  earth  model  as  a  series  of 
radially  inhomogeneous  layers,  each  having  a  constant  raidal  gradient  in 
velocity  and  density.  Seismogram  synthesis  can  then  be  most  efficiently 
achieved  using  the  propagator  matrix  solution  in  either  a  reflectivity  or 
a  mode  summation  method. 

Solutions  to  Potential  Wave  Equations. 

The  Fourier-transformed  potential  wave  equations  satisfied  by  P,  SV, 
and  SH  scalar  potentials  are  given  by  Richards  (1974)  as 
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0  2  K,  T.2 

V2P  +  ■  -  P  +  e  P  =  — i- 
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2  K 

y  ,  V  =  -  _ 
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V2H  +  -2^-  H  +  e  _  H  =  0 
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1  f  32p  8  ,BV.l 

2  L^2'3'  r  J 

\  fap  bv1 

„2  L3r  "  r  J 


(  -la) 


(  -lb) 


(  .lc) 


for  P,  SV,  and  SH  waves  respectively,  where  B  is  an  operator  acting  as  the 
non-radial  part  of  the  Laplacian.  Separation  of  variables  allows  the 
solutions  to  be  expressed  as  a  sum  over  Legendre  polynomials: 


P(r,u)  =  ^  ]  w(r,  n)  P^cosA) 


(  .2a) 


V(r,w)  =  ^  J  [x(r ,n) / (iwp) ]  Pn(cosA)  (  .2b) 


H(r ,w) 


=  y(r,n)  pn 


(cosA) 


(  -2c) 


The  constant  factor  iWp  (i  =  /-[  ,  w  is  radial  frequency,  p  ray 

parameter)  divides  the  radial  eigenfunction  for  SV  waves  to  recognize  its 
dimensional  difference  from  that  of  the  radial  eignefunction  for  P  waves. 
This  difference  arises  from  the  definition  of  P  and  SV  potentials  and 
must  be  considered  when  determining  the  order  in  frequency  of  their  coupling 
terms. 


L 


Substituting  (  .2a-  .2c)  in  (  .la-  .lc)  results  in  equations  for  the 


radial  functions: 
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(  -3c) 


where  W  ,  X  ,  Y  equal  rw  ,  rx/ (imp)  ,  ry  respectively. 

Taking  a  Liouville  transform  of  eqs.  (  .3a-  .3c),  making  the  variable 
change  proposed  by  Langer  (1949),  and  substituting  up-1/2  for  n  gives 
the  equations: 
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The  variables  £  and  functions  0  ,  0  are  defined  by 

* » S  r  s 

(  .4f) 

(  .4g) 

where  r^  and  rg  are  the  respective  turning  point  radii  for  P  and  S  waves 
in  the  medium.  The  functions  and  are  the  Schwarzian  derivatives 

resulting  from  the  Liouville  transformation  plus  a  term  reuslting  from  the 
substitution  for  n  : 


r  _  9”’  +  3.0H  +  -A_ 

5  -  201-  +  4(0'  )  +  .  2 

4r 


The  Schwarzian  derivative  can  be  evaluated  in  terms  of  the  physical  variables 
r  ,  a  ,  8  ,  and  p  by  the  identity 
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where  the  subscripts  P  and  S  are  omitted  for  evaluation  in  either  the  P 
or  S  velocity  profiles  and 


T  = 


dr 


The  radial  eigenfunctions  may  now  be  found  to  any  desired  order  in 
frequency  using  the  fundamental  series  solution  given  by  Olver  (1976) . 
For 
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where  Ai  represents  any  Airy  function  solution  of  the  zeroth  under 

equation  and  A  ,  B  ,  A  ,  B  are  functions  to  be  determined  by 
m  m  m  m 

substitution  of  the  series  of  eq.  .4a  and  ,4b  and  equating  terms  of  equal 
order  in  frequency.  Thus  it  can  be  determined  that  Aq  is  always  constant 
in  radius  for  P,  SV,  and  SH  waves  and  in  terms  of  physical  variables  Bq 
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where  opS  ,  0gp  are  coupling  coefficients  defined  by 
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in  which  the  subscript  (o)  denotes  the  zeroth  order  term  in  the  solution 

series  for  W  and  X  .  After  substitution  of  W  and  X  above,  it 

o  o 

can  be  shown  that  the  coupling  coefficients  Opg  ,  Ogp  have  a  strength 
proportional  to  p/r  .  Thus  at  vertical  incidence  on  a  region  of  intense 
vertical  gradient  the  coupling  coefficients  0^  ,  0gp  -  0  .  This 

continuum  property  mimics  the  non-conversion  of  P  and  SV  waves  vertically 

. . . 


incident  on  a  discontinuity.  For  non-vertical  incidence  the  coupling 
coefficients  introduce  a  phase  factor  determined  by  the  strength  of  the 
function  K^/p  between  r  and  the  turning  point  radius.  The  phase 
factor  acts  to  correct  for  the  phase  accumulated  along  the  path  of  a  con¬ 
verted  wave  type  from  the  conversion  point.  Because  the  coupling  coefficients 
require  the  evaluation  of  additional  Airy  functions  in  the  integrand  of 
integral  for  the  Bq  term,  the  efficiency  of  the  higher  order  potential 
solution  in  seismogram  synthesis  would  be  poor. 

Propagator  Solution  for  Wave  Propagation. 

An  alternative  description  of  seismic  wave  propagation  in  an  inhomo¬ 
geneous  medium  can  be  made  in  terms  of  the  matrix  equation  satisfied  by 
components  of  displacement  and  stress  (Gilbert. and  Backus,  1966).  The 
wave  equation  is  thus  written  in  the  form 

3F 

~  =  sM(r)  FJr) 

where  is  a  funcamental  matrix  solution  whose  elements  are  Fourier- 

transformed  components  of  displacement  and  stress.  (s  =  -iw  for  a  forward 
Fourier-transform  sign  of  the  form  /e+*Wt  dt) . 

In  obtaining  higher  order  solutions  for  F  Chapman  (1973)  decomposed 
the  matrix  M  into  a  sum  of  matrices  of  differing  order  in  frequency,  i.e.. 


J 


j=0 


where  J  =  2  for  SH  waves  and  4  for  P  or  SV  waves.  Wasow  (1965)  described 
how  a  uniformly  asymptotic  solution  to  such  a  system  can  be  formulated 
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1 


similar  to  the  solution  eqs.(  .5a-  .5b)  derived  by  Olver  (1976)  for  a  single 
component  system.  Using  a  series  of  similarity  transformations  and  hanger's 
(1949)  variable  change,  Chapman  (1974)  obtained  matrix  equations  analogous 
to  the  transformed  equations  (  .4a-  .4b), 


/v 


(  .8) 


and  solution  formulae  analogous  to  eqs.  (  .5a-  .5b), 


(  -9) 


where 


rAi(£S2/3) 


Bi(£S2/3) 


S~2/3  Al(£S2/3)  S_1/3  Bi(£S2/3)i 


(Here  the  symbol  'v  over  a  matrix  denotes  that  the  Langer  variable  change 
of  eqs.  (  .4f-  .  4g)  has  been  made.)  Chapman  (1974)  determined  the  matrix 
in  a  flattened  earth  as 


(  .10) 


where  I  is  the  identity  matrix  and  is  a  reference  depth  most  conveniently 

chosen  as  the  earth's  surface. 
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(K) 

Using  Chapman's  (1974)  results  for  higher  order  matrices  Mv  and 
Wasov’s  (1965)  theorems,  the  matrices  etc.  can  be 

determined.  The  procedure  solves  for  the  functions  »  9^^  *n 

solutions  of  the  form 


,(0) 


=  q 


(0) 


(£)  I  +  q 


(0) 


(S)  h(0) 


(  .lla) 


,(1) 


+  9™ 


a) 


i(1) 

l2 


a)  m(0) 


c  .nb) 


subject  to  the  compatibility  conditions 


M(0)  M(0)  =  0  for  K  =  0 

K 

=  D(K-1) '  -  £  M(j)  D(K'j)  for  K  *  0  (  .12) 


j=l 


where  is  a  particular  solution  and  ( ’)  signifies  differention  w.r.t. 

_P 

the  Langer  variable  £  .  For  example  by  considering  the  matrix  system 

for  SH  waves  a  particular  solution  to  the  compatbility  condition  for 
is  given  in  a  flattened  earth  by 


(  .13) 


where  (•)  signifies  differention  w.r.t.  to  depth  Z  and  C  is  a  constant. 

By  again  substituting  the  general  solution  (lib)  for  in  the  compatability 

condition  it  can  be  determined  that 
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(  .14a) 


and 


12 


+  m 


22 


d21 

v/iT 


21 


dZ 


(  .14b) 


where  a  is  the  radius  of  the  earth  and  m^  ,  d^.  denote  the 

ij  elements  of  the  and  Dp^  matrices  respectively.  Now  substituting 

this  result  in  eq.  (  .9)  with  expressions  for  the  elements  of  , 

the  solution  for  the  fundamental  matric  F  to  second  order  becomes 


2o/s 


(  .15) 


where 


When  the  quantity  d^2  =  q^  is  expressed  in  a  single  integral  IgH 
over  depth,  the  quantity  d2^  +  £  q^^  expressed  in  terms  of  ISH  ,  and 
the  varible  of  integration  changed  from  the  flattened  depth  coordinate 
Z  to  radius  r  it  can  be  shown  that 


£T  +  6S  dr  _  ~ISH 

V  2Qs 


(  .16a) 
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and 


21 


+  5q 


(1) 


0  T  •  •  • 

-1  ^  S  sH  1  JJ  +  i 

|3/2  2  2p  • 


(  .16b) 


The  results  given  by  (  .16a-  .16b)  for  a  spherical  layer  are  equivalent 

to  those  obtained  by  Woodhouse  (1977)  for  SH  waves  in  plane  layer.  Note 

that  the  integral  IgH  is  the  same  as  that  derived  for  the  higher  order 

potential  term  Bq  in  eq.  .6c  for  SH  waves. 

Practical  Evaluation  of  Seismic  Displacement 

By  Watson  transforming  the  partial  wave  expansion  the  seismic  dis- 
P  SV  SV 

placement  components  u  ,  u  ,  u  can  be  represented  in  the  frequency 
domain  as 


io,3 


ur  (r,  &o>  <J>,w)  = 


4frp 


M°  />  ■ 
T  /  F  w 
a 


(1)(rQ)  w’(1)(r)  Q(2)  dP 


(  .17a) 


uf  (ro’Ao’  U) 


iw3  M 


‘•"Wo  r 


f  FSV  xU)(ro)  x'(1)(ro)  Q(2>dp 


(  -17b) 


ufH  (r’  A0’  *■  w)  =  f  fSH  y(1)(0  y’(1)<0  Q(2)  dP  (  -i7c) 

“  °  /.  trr\  O*  J 


4Tfp6 

O  O  i 


for  a  shear  dislocation  centered  at  radius  r  with  seismic  moment  M 

o  o 

P  SV  SH 

and  a  radiation  pattern  described  by  the  factors  F  ,  F  ,  and  F 
For  a  continuous  earth  with  regions  of  intense  gradient  the  last  semi-annual 
report  considered  the  feasibility  of  calculating  displacement  by  inte¬ 

gration  in  the  complex  p  plane  using  the  respresentations  given  by 
eqs.  (  16a-  .16b).  If  the  regions  of  intense  gradient  are  sufficiently 


deep,  the  Airy  functions  in  the  higher  order  representation  for  the  radial 
eigenfunctions  w,  x,  y  may  be  replaced  by  their  asymptotic  approximations. 
Thus  w,  x,  y  are  approximated  by 


(1)  1/2  +i7r/4  +T 

<2)_.ao  6  *"  P 

7  ^  - - - 

rwQ 


1  +  —  [— - - — ] 

-  u  1  2  72x  J 


(2)  g1/2  e+i^  e±iTS 

x  '  o 


-iojp 


rui  Q. 


i  +  A  -  -L  ] 

-  u  1  2  72tJ 


(1)  1/2  +i»/4  +ix 

P e  e—  o 


y<2>«^ 


rco  Qc 


1  [I 


-  w  1  SH  72  r0 


(  .18a) 


(  .18b) 

(  .18c) 


for  (1)  up  or  (2)  downgoing  waves. 

The  representations  given  by  eqs.  (  .17)  and  (  .18),  however,  are 
additionally  limited  in  practical  use  to  situations  in  which  the  coupling 
coefficient  present  in  ,  I_„  ,  and  I  „  can  be  neglected.  The 

r  on  bv 

coupling  coefficients  themselves  involve  the  evaluation  of  Airy  functions 
of  both  small  and  large  arguments,  making  the  evaluation  of  Ip  and  Igy 
very  inefficient.  A  far  better  representation  in  practical  use  for  P-SV 
displacements  is  given  by  the  higher  order  fundamental  matrices  determined 
by  Woodhouse  (1977) .  The  higher  order  matrix  correction  to  the  fundamental 
matrix  for  P  and  SV  waves  fully  accounts  for  wave-type  conversion  and  in¬ 
volves  integrals  no  more  complicated  than  that  of  the  I  type.  A  given 
displacement  component  observed  at  the  earth's  surface  can  be  calculated 
from  the  product  of  propagator  matrices  for  each  model  layer  with  velocity 
and  density  describable  by  functions  analytic  in  velocity.  For  example, 
the  propagator  matrix  from  the  boundary  radii  r ^  to  r^  of  such  a  layer 


can  be  written  as 
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P(r2,  rx)  =  F(r2)  F-1^)  (  .19) 

%  r"1  R(r2)  (1  +  iuf1  r2  T(1)(r2))  |(r2>  G(r2) 

*  §~1(r2)  ri  I(1>(r1))  i’1^)  T1 

Woodhouse  (1977).  The  elements  of  R  are  simple  functions  of  elastic 
constants  density  and  ray  parameter.  For  practical  evaluation,  the  matrix 
products  <£  G  and  G  ^  <J>  ^  can  be  expressed  as  matrices  containing 

up  and  downgoing  generalized  cosines  and  radial  eigenfunctions  using  the 
definitions  of  Richards  (1976) 


where 


cosj 


+  iwy 


(2) 
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Fuchs  and  Muller  (1971)  describe  how  a  displacement  component  observed 
at  the  earth’s  surface  can  be  calculated  using  an  integral  representation 
equivalent  to  that  eqs.  (  .17a-  .17c)  but  in  which  a  reflectivity  function 
calculated  from  a  product  of  propagator  matrices  substitutes  for  the  radial 
eigenfunctions  w  ,  x  or  y  .  An  alternative  representation  of 
displacement  does  not  Watson  transform  the  partial  wave  series  and  instead 
takes  a  truncated  sum  of  the  discrete  modes  n  ■*  0)p  -  1/2  (Sato  and 
Usami,  1963).  This  representation  allows  synthesis  of  a  longer  portion  in 
time  of  the  seismogram  from  body  to  surface  waves.  The  use  of  Woodhouse's 
(1977)  propagator  matrices  in  either  method  eliminates  the  need  for  describing 
the  earth  model  as  a  stack  of  many  plane  homogeneous  layers.  Langer's 
approximation  to  radial  eigenfunctions,  embedded  in  the  Airy  functions 
in  the  <j>  G  and  G  ^  <j>  *  matrices,  fully  accounts  for  velocity  gradients 

and  boundary  curvature  of  radially  inhomogeneous  layers.  A  matrix  solution 
method  is  capable  of  handling  mode  conversion  more  efficiently  than  the 
potential  solution. 

Since  layers  may  be  inhomogeneous,  a  question  then  remains  of  how 
many  inhomogeneous  layers  are  necessary  to  describe  an  earth  model.  In 
principle  the  earth  model  may  be  completely  continuous.  At  sufficiently  high 
frequencies  the  next  higher  order  term  for  the  fundamental  matrix  accounts  for 
the  effect  of  regions  in  which  velocity  and  density  gradients  are  large 
and  rapidly  changing.  In  practice,  however,  it  is  seen  that  it  is  best  not 
to  represent  a  region  of  rapidly  changing  velocity  and  density  gradient  as 
a  single  inhomogeneous  layer,  e.g.,  representing  a  low  velocity  zone  by 
a  parabolic  function  or  a  transition  zone  by  a  hyperbolic  tangent  function. 

In  these  cases  the  integrand  of  the  integrals  I  etc.  in  the  higher  order 

art 
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matrix  would  be  large  and  rapidly  changing  at  those  depths  at  which  the 

velocity  and  density  gradients  are  rapidly  changing  (Figure  8) .  One  must 

choose  between  either  increasing  the  number  of  inhomogeneous  layers  in  the 

earth  model,  and  thereby  the  number  of  fundamental  matrices  to  be  evaluated, 

or  increasing  the  computation  time  and  magnitude  of  the  higher  order  correction 

to  the  fundamental  matrices.  Representation  of  the  earth  model  as  inhomogeneous 

layers  of  nearly  constant  depth  gradient  in  velocity  and  density  clearly 

most  efficiently  compromises  between  these  choices.  This  representation 

minimizes  the  effect  of  the  higher  order  correction  to  the  fundamental 

matrix  and  in  many  cases  allows  the  integrals  I  etc.  to  be  evaluated 

SH 

using  the  three  point  evaluation  scheme  described  by  Jeffreys  and  Jeffreys 
(1956). 

In  the  next  research  period  a  practical  computation  code  will  be  developed 
Incorporating  higher  propagator  matrices  in  a  either  a  reflectivity  or  mode 
summation  method  of  seismogram  synthesis.  The  method  will  be  tested  in  speed 
and  accuracy  against  other  exis^  methods  of  seismogram  synthesis. 
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Appendix  2:  General  Elastodynamics  Representation  Theory  for  Inhomogeneous 
Media  with  Moving  and  Fixed  Internal  Boundaries. 

(1)  Conservation  Equations  in  a  Linear  Elastic  Medium 

The  linearized  equations  of  motion  may  be  written  as+  (Archambeau  and 
Minster,  1977) : 

l>aY»»Y  -  pfo  ;  ot.Y.fi  -  (1,2, 3, 4)  (1) 

Here,  u^  is  the  space- like  displacement  field 

UY  B  (VW0) 

while  f  is  the  space-like  force  density 

fa  = (fi»f2>f3»°) 
and  x  is  the  four  vector: 

Xy  ~  =  ^xi»x2*x3,t^ 

The  elastic  operator  is  given  by: 


+The  summation  convention  for  repeated  indices  applies  throughout..  Indices 
excluded  from  this  rule  will  be  enclosed  in  parentheses,  e.g.  ,  G^j  Uj  implies 
no  sum  on  n,  but  summation  on  the  j  index  over  its  range  of  values  (1,2,3). 

Also,  Latin  indices  will  range  from  1  to  3,  Greek  indices  from  1  to  4. 

Cartesian  tensors  are  used  throughout,  so  while  subscripts  and  superscripts 
are  used  as  is  notationally  convenient,  they  do  not  denote  covarient  or 
contravarient  components . 
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with  CagY^  the  elastic-inertial  tensor: 


“a0Y«* 


~  ~  (1*2,3) 


-a0Y«  "  uijkl 


Ci4k4  °  C4i4k  “  p6ik;i,k  “  (1,2,3) 


Co8Y«  “  °  J  otherwise 


(3) 


with  the  usual  elastic  tensor.  For  an  isotropic  elastic  medium 


Cijkl  =  A6ij6lk  +  w(6il6ik  +  6ik6il)* 


^ijkl  an<*  Cq0y^  °^ey  t^ie  same  symmetry  conditions  in  all  cases,  namely: 


Cijkl  =  Cjikl  =  Cijlk  Cklij 


and 


CoiBy^  “  ^BoyS  =  ^oB6y  =  S'ScxB 

An  alternate  form  of  the  equations  of  motion  is  obtained  by  defining 
the  generalized  inertial-stress  tensor: 


9uy 

TaB  ~  ^oBy^ 


(4) 


rr 


Then,  from  (1),  we  have  the  equation  of  motlon+ 


Ta$,8  “  pfa 


Finally,  the  natural  boundary  conditions  in  elastodynamics  may  be 

expressed  in  terms  of  continuity  conditions  involving  r  Thus,  if  we 

op 

consider  general  media  in  which  we  admit  the  possibility  of  moving 
boundaries,  such  as  a  growing  failure  zone  or  moving  phase  boundaries,  etc., 
then  the  boundary  conditions  can  be  expressed  in  compact  form  if  we  define 
a  "space-time  normal"  to  the  surface  as  a  four  vector  n0,  where 

P 


hg  —  ,n2  »n^ j  —  UJ  n^) 


where  n  =  (n^.n^.n^)  is  the  ordinary  spatial  normal  to  the  boundary  surface 
and  U*  is  defined  as 


u*  -  u*  "  v* 


where  v  is  the  particle  velocity  within  the  medium  into  which  the  normal 
to  the  surface  is  directed  and  U  is  the  velocity  of  the  boundary  surface. 
The  boundaries  within  and  enclosing  the  medium  will  be  designated  by 


+The  usual  notational  convention  for  partial  derivatives  will  also  be  used 
occasionally,  that  is: 


8eTa0  =To8>B  "  3x„ 


I 
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the  symbol  E.  The  boundary  condition  expressing  conservation  of  momentum 
across  any  interface  or  exterior  boundary  of  the  medium  is  (Archambeau 
and  Minster,  1977)+ 

-  c  (8) 

written  out  in  component  form  this  becomes,  using  the  definitions  of  t 

ap 

and  n^: 

-  V  “A  -  0  (8-a) 

where 

vi  "  v£  “  U£. 

is  the  particle  velocity  relative  to  the  moving  boundary.  The  tensor 

Tjj^  is  the  ordinary  Cauchy  stress  tensor. 

Clearly,  if  the  boundary  moves  with  the  particles  of  the  medium,  in  the 
*  * 

sense  that  U^n^  =0,  then  v^  n^  »  0  also, and  the  equations  (8)  reduce  to 
the  regular  continuity  of  traction  conditions  across  E,  that  is: 

nv,nr  - 0  <*-« 

The  double  bracket  notation  IMe  is  used  to  denote  the  change  in  a 
function  F  on  crossing  a  surface  E.  That  is: 

HFlIj.  -  Ftfj)  -  F(E2) 
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The  remaining  boundary  conditions  express  conservation  of  mass  and 
energy.  They  are  (Archambeau  and  Minster,  1977): 

C  ~  0  (9) 

and 

□XpEv*  ~  vkTki  +  <!*)  ni^2:  “  0  <10> 

Here  q  is  the  heat  flux  and  E  is  the  total  energy 

pE  =  pu  +  p/2v^v^  +  p«J> 


where  u  is  the  specific  internal  energy  and  <f>  the  body  force  density  potential, 
3 


so  ,f  = 
a 


3x 


If  the  boundary  moves  with  the  material  so  that  the  boundary  is  carried 

'fc 

normal  to  itself  with  the  particles  (i.e. ,  U^n^  =0),  then  (9)  reduces  to 
continuity  of  the  normal  particle  velocity  across  E,  so: 


=  ° 


(9-a) 


In  this  case  (10)  reduces  to 


(10-a) 


where  t^  ■  ^kjj,n£  are  t^ie  comP°nents  of  traction.  This  relation  simply 
expresses  the  rate  of  heat  production  on  the  boundary  (in  terms  of  a  jump 
in  heat  flux)  if  there  is  relative  slip  of  the  material  on  opposite  sides 
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of  the  boundary.  This  slip  is  constrained  to  be  along  the  boundary, 

/ 

however,  because  of  the  constraint  imposed  by  (9-a). 

If  the  physical  boundary  is  such  that  a  no  slip  condition  is  warranted  - 
that  is  lTv^Tl^  =>  0  -  then  (10-a)  reduces  to: 

□Vplj.  -  0  (10-b) 

This  just  states  that  the  heat  flux  must  be  continuous  across  the  boundary. 

Since  thermal  effects  are  of  second  order  in  such  a  situation  ('welded" 
boundaries) ,  then  this  boundary  condition  may  be  neglected  in  solving  for 
the  field  u^  from  the  equations  of  motion.  The  boundary  condition  of 
importance,  in  addition  to  (8-b) ,  is  then  the  condition  of  no  slip,  which 
requires : 

iivknr  *=  °  ci!) 

If  this  continuity  condition  is  satisfied, then  (9-a)  is  automatically 
satisfied.  Usually  (11)  is  expressed  more  strongly  through  the  use  of  the 
Sufficient) displacement  continuity  condition.  In  particular 

(Xu^TIg  -  0  (11- a) 

It  is  important  to  note  that  (11)  and  (11-a) ,  amount  to  assumptions 

regarding  the  physical  processes  at  the  boundary  and  are  not  necessary 

conditions,  the  most  general  case  being  covered  by  the  relations  (8),  (9)  and  (10). 

However,  it  is  often  the  case  that  the  boundary  may  be  considered  to  be 
"welded"  (for  solid-solid  interfaces)  and  then  (8-b)  and  (11)  are  the 
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appropriate  boundary  conditions. 

In  the  cases  to  be  treated  here  we  are  concerned  with  the  simpler 
case  of  welded  internal  boundaries  (solid-solid  contacts).  In  this  case 
we  consider  solutions  of  the  system: 


0[8  u  =  0 
ay  y  l 

uu^rij. » o 


(12) 


where  the 


"boundary  operator" 


is  defined  as: 


(13) 


On  external  boundaries  or  solid-fluid  interfaces  the  displacement  condition 
is  either  absent  or  replaced  by: 


(14) 


4  JL 

if  one  or  both  materials  are  fluid.  Since  the  condition  that  U„n„  =  -  v,n.  **  0 

Z  Z  Z  Z 

applies  (the  motion  of  the  boundary  normal  to  itself  is  with  the  material 
particles)  then  the  boundary  conditions  in  (12)  reduce  to  the  ordinary 
continuity  conditions  for  traction  and  displacement. 


•J' 

In  the  linearized  elastodynamic  theory , boundary  motion  is  neglected  in  the 
sense  that  all  fields  are  evaluated  at  the  undeformed  boundary  position. 

For  an  external  boundary  this  leads  to  neglecting  the  continuity  of  normal 
velocity  condition. 
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The  boundary  conditions  expressed  in  (12)  may  also  be  written  in 

/ 

the  matrix  operator  form: 


(; )  -1 


(2)  Green’s  Tensors  for  the  Linear  Elastic  Medium 

The  Green’s  tensor  associated  with  the  elastic  operator  L  can  be 

ay 

defined  by: 


L  G®(x;x  )  =  A^(x,x  ) 
ay  y  —  — o  a  —  — o 


Here  A^  is  the  generalized  delta  function 


Ao  H  4ir 


(16~a) 


As  before  the  coordinate  vector  x  is  the  four  vector  with  components 

8 

(x. ,  x_,x _,t).  Similarly,  x  is  a  coordinate  four  vector.  Thus  G  is  a 
x  4  ^  “O  Y 

g 

second  order,  two  point  tensor.  Specifically,  G^  corresponds  to  the 
space-like  displacement  at  x  in  the  6  direction  due  to  a  point  (vector) 

g 

force  at  3^  in  the  y  direction.  Since  G  is  space-like,  then  the  time-like 

4  3 

component  G^  is  identically  zero.  G^  plays  the  roles  of  a  propagator 
for  elastic  displacement  fields  since  solutions  of  (15)  produce  functions 

g 

G“  obeying  the  relations  required  for  displacement  field  propagation  in  the 
medium  from  one  point  (x^)  to  another  (x).  In  view  of  the  meaning  to  be 

g 

attached  to  the  Green’s  function  Gy,  only  causal  Green's  functions  are 
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desired,  so  that  G  must  also  satisfy 
Y 


A  4  4 

G^(x,x  )  *•  0  ;  x  <  x!: 

y - o  u 

8 

This  means  that  G^  must  have  the  property 


G^(xk,x*;x£,XQ)  =  G®(xk,-XQ;xk,-x^) 


The  transmission  properties  of  the  medium  are  expressed  by  the  opera' 
However,  the  full  specification  of  the  medium,  including  the  nature  of  the 
external  and  internal  boundaries,  is  achieved  by  specifying  boundary 
conditions. 

Therefore  complete  specification  of  requires  a  statement  of  boundary 
conditions.  For  compatibility  with  the  boundary  conditions  of  (12) 
we  can  take: 


aV\  nE  -  o 

cig®d£  -  o 


or 


(17) 


(17-a) 


At  "unwelded"  boundaries  (solid-fluid;  fluid,  solid-vacuum)  only  the  first 
of  the  relations  in  (17)  applies  in  general,  while  the  second  condition 
is  replaced  by  the  normal  velocity  continuity  condition. 


89 


“°;.*'VDt  -  o 


(17-b) 


when  a  fluid  is  involved.  At  the  exterior  boundaries  of  the  medium  (usually 
considered  as  an  interface  with  a  vacuum)  only  the  traction  condition  need 
be  applied  (the  first  condition  in  (17)). 

The  equations  (15)  and  boundary  conditions  like  (17)  plus  the  causality 

g 

condition,  define  the  complete  operator  for  G^. 

The  operator  to  be  used  to  generate  the  Green’s  tensor  ^(xjx^) 

associated  with  the  displacement  field  u„(x)  is  ideally  to  be  chosen  so 

P  — 

g 

that  G^  can  be  employed  to  propagate  known  initial  values  or  boundary 

values  of  the  field  u g  through  the  medium,  so  that  Ug  can  be  predicted  at 

g 

other  space-time  points.  Thus,  G^  is  to  be  such  that  it  acts  like  a  simple 
transfer  function.  A  necessary  condition  for  G^  to  act  in  this  manner  is 
that  it  satisfy  (15),  plus  the  causality  condition,  plus  boundary  conditions 
accurately  reflecting  the  properties  of  the  medium  like  those  of  (17).  However, 
as  will  be  seen  in  the  following  section,  when  the  field  ug  is  known  on  a 
particular  boundary  then  it  is  appropriate  to  use  homogeneous  boundary 

g 

conditions  for  G^  on  this  boundary.  (This  statement  applies  to  time-like 
boundaries,  for  initial  values,  as  well  as  to  space-like  boundaries.) 
Consequently,  the  space-time  boundary  condition  for  Eq,  the  particular 
surface  over  which  the  value  of  u^,  is  known  or  specified,  is  of  the 
homogeneous  form, 

B  G8  L  =0, 

Y 
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with  the  notation  here  specifying  that  the  quantity  is  evaluated  as  a 

limit  from  within  the  medium  on  one  side  of  the  surface  1  .  The  role 

o 

g 

of  this  choice  of  boundary  condition  on  for  will  be  apparent  in 

g 

the  following  section,  where  it  will  be  seen  that  will  then  have  the 
characteristics  of  a  transfer  function  for  the  displacement  field  ue. 

p 

However,  it  is  both  necessary  and  advantageous  to  use  a  causal  Green’s 
function,  from  the  wide  class  generated  by  (15) ,  satisfying  boundary 
conditions  other  than  those  of  the  "associated"f ield  uQ.  This  arises  from 

P 

the  fact  that  it  is  usually  as  difficult  to  determine  the  appropriate 

Green's  function  (satisfying  all  the  boundary  conditions)  as  it  is  to 

obtain  complete  solutions  for  u.  itself.  Thus  it  is  usual  to  relax  the 

p 

boundary  conditions  on  the  Green's  function  and  to  obtain  a  Green's 

function  that  provides  an  approximate  "transfer  function"  for  u0.  Such 

p 

procedures  result  in  approximations  for  the  field  uD  elsewhere  in  the 

p 

space.  Aspects  of  this  approximate  technique  are  also  discussed  by 
Archambeau  and  Minster  (1977),  for  elastodynamic  relaxation  source 
representations . 


The  formally  defined  properties  of  the  Green's  function,  in  particular 


the  differential  relation  (15) ,  can  be  used  along  with  the  differential 
equation  involving  u„  in  (12)  to  obtain  an  integral  equation  relating  the 

p 

value  of  the  field  u  at  all  points  in  the  space  to  known  or  specified  values 

P 

on  space- time  boundaries. 

In  particular,  using  Green's  theorem  in  the  four  dimensional  space  ft, 
we  have  for  two  fields  with  appropriate  properties  of  continuity  (Archambeau 
and  Minster,  1977): 

(4t,v)ft  *  (w»l*v)R  +  /J6>0d4x  (18) 

where  the  inner  products  are  defined  according  to: 


(Lw,v)n  =  fv  L  w  d^x  (13-a) 

—’—  ft  *  a  an 

* 

Here  (L^)  =  L  is  defined  in  (2).  The  operator  L  is  the  adjoint  operator 
to  L.  It  can  be  easily  shown  that  1  is  a  self-adjoint  differential  expression 
that  is:  L  =  L  . 

The  quantity  Jft  is: 


va  Ca8Y5  wy,6  Wy  CaSy6  Vy,6 


(18-b) 


with  C  „  .  the  elastic-inertial  tensor  defined  earlier.  It  is  evident  that 
a8yo 

a  formal  application  of  Gauss'  theorem  to  the  final  integral  in  (18)  produces 
a  "surface"  integral  over  the  boundary  of  ft,  involving  the  projection  of  J- 

P 

on  the  normal  to  this  surface  as  the  integrand.  We  note,  in  fact,  that  this 
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is  J-,n,, ,  and  that 

P  P  / 


Vs 


v  [B  w  1  -  w 
a  ay  y  a 


[B  v  ] 
«Y  Y 


where  B^JLs  the  boundary  operator  of  (12)  and  (17).  It  is  this  relation  that 
explicitly  displays  the  role  of  the  boundary  conditions  in  the  integral 
equation  (or  solution)  for  u^,  which  will  emerge  from  (18). 

In  the  important  cases  in  which  the  fields,  or  L  itself,  have  first  or 
higher  order  discontinuities,  then  (18)  cannot  be  employed  in  a  direct, 
straightforward,  fashion.  The  physical  situation  in  which  the  medium  has 
internal  boundaries  across  which  the  material  properties  change  is  a  particularly 
important  example.  Further,  this  discontinuous  behavior  may  involve  either 
fixed  or  rapidly  moving  boundaries,  so  that  it  may  be  space-time  coupled. 

Thus  it  is  critical  for  applications  that  the  more  general  discontinuous 
case  be  considered  explicitly.  To  do  so  we  can  redefine  the  inner  product 
over  the  space  ft,  as: 


(L»,v)a  -  f 

Jm 


v  L  w 

a  ay  Y 


d*x 


+  . . .+ 


/ 

{m 


v  L  w 
a  ay  Y 


d4x 


(19) 


(P) 

where  the  ft  are  sub-regions,  divided  along  space-time  boundaries,  within 
which  all  quantities  in  the  integration  are  continuous.  Implied  here  is  a 
separation  of  the  sub-regions  by  a  space-time  strip  of  width  n2e"  along  the 
surface  of  discontinuity,  with  the  limit  e -*■  0  applied  to  the  sum  of  integrals. 


We  can,  however,  also  write  this  as 
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(Lw,v) 


f  4  f 

n=l  v  L  w  d  x  =  /  v  L  w 

ft  I  a  ay  Y  I  a  ay 

:{i).  „(n) 

ft  ...  7  i 


d*x 


(20) 


if  we  merely  take  note  of  the  fact  that  this  integral  form  can  always  be 

written  as  a  sum  of  integrals  over  the  disjoint  sub-regions  defined  by 

fiOEj  (or  ft^^©-*©fl^),  where  E  are  the  internal  surfaces  of  discontinuity 

within  0  ,+  To  be  definite  ft^  is  bounded  by  one  of  the  exterior  boundaries 
(N) 

of  ft  while  ft  has,  as  a  boundary  any  other  exterior  surface.  (This  "surface" 
may  be  at  the  origin  or  at  infinity.)  Then  all  the  ftv  ,2  4P$H,  are 
bounded  by  Internal  surfaces  of  discontinuity  in  the  medium. 


With  this  inner  product  definition  in  ft  we  can  apply  the  generalized 

(py 

Green's  theorem  of  (18)  to  each  subregion  ftv  .  It  follows  from  (19)  and 
(18)  therefore,  that: 

'•  (w,L*v)fl  +  J ?  dSc  (21) 

001 

Here  the  inner  products  are  defined  by  (19)  or  (20)  and  the  previously 
mentioned  limit  procedure  applied  to  the  subregion  integrations  is  implied. 


+The  symbols  ©  and  9  denote  the  set  theoretic  sum  and  difference.  Thus, 
ftQlj  denotes  the  space  ft  with  the  points  of  all  the  surface*  Ej-  deleted. 
For  N  such  surfaces,  this  is  equivalent  to  ft  ©••©ftW. 
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This  result  is  nearly  the  same  as  (18) ,  and  indeed  is  the  same  if  we 

/ 

always  think  of  ft  as  the  set  of  subregions  in  which  the  fields  and  the 
differential  operator  are  continuous.  However,  (21)  makes  this  fact  clear 
and,  in  addition  to  being  more  explicit  than  (18)  in  showing  precisely 
how  internal  boundaries  of  discontinuity  enter  the  problem,  it  is  also 
rigorously  derived  for  the  general  case. 

Some  care  must  clearly  be  used  in  applying  Green's  theorem  along  with 
Gauss'  theorem  to  the  elastodynamic  problem  to  obtain  an  integral  relation 
for  the  displacement  field  in  view  of  the  fact  that  internal  boundaries 
are  important.  To  obtain  the  desired  relation  we  take  w^  in  (21)  to  be  a 

O 

two  point  causal  Green’s  tensor  Gp  satisfying  (15)  and  v  to  be  a  displacement 

Y  ® 

field  u^  satisfying  the  differential  expression  in  (12).  Then  considering  the 

inner  product  relation  (21)  over  the  coordinate  space  (the  "source 

coordinates,  instead  of  the  "observer"  coordinates  x.)  along  with  differential 

expressions  with  x  as  the  independent  variable,  one  has  (Archambeau  and 
o 

Minster,  1977)+: 


+Archambeau  and  Minster  denote  the  integration  region  simply  as  ft,  and 
implicitly  use  the  general  definition  of  ft  as  the  set  of  subregions  in 
which  continuity  holds.  Here,  however,  we  will  be  more  explicit  about 
this  restriction  on  the  integration  and  write  ftSE^,  with  the  implied  limit 
to  be  taken. 
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4nUy(x)  =  J 


p<*o>  <£(x;x  )  * 


SIQZ, 


(22) 


-/b  l 

rising  l 


GP(x;x  )tq(x)-u(x)  GP0(x 
a  — ’— o  ag  — o  a— o  ag  — 


[x;x  )  1  d*x° 

--x,  j 


The  displacement  field  term  on  the  left  side  arises  from  an  integration  of  the 
left  side  of  equation  (21)  using  the  generalized  delta  function  in  the 
differential  equation  for  the  Green's  function.  The  first  tern  on  the  right 
side  is  the  same  as  the  inner  product  appearing  in  (21),  with  1**^^  ■  L^u^ 
replaced  by  its  equivalent  pfQ,  from  equation  (12).  Similarly,  the  final 
integral  is  as  given  by  (21) ,  with  ^  being  replaced  by  its  substitutional 
equivalent,  D,  where 

p  ,p 


Jr(G,u) 


■Gw(x;x  )td(x)-u(x)  GP(x;x  ) 
a  — *-o  ag  -o  a'— o'  ctg  —  ~o 


(23) 


The  Green's  inertial  stress  tensor  Go  is  defined  analogously  to  the  regular 

ap 

inertial-stress  tensor  (equation  4);  in  particular: 


g!!d(x;x  )  =  c 


ag'-’^o 


agy<5 


(^>  3Z - 


(24) 


As  with  JgOg,  we  note  that: 


TaB*S 


B  u  ’V 

ay  Y  I 

B  G*  ) 
ay  y 


(25) 
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so  that  the  previously  discussed  boundary  conditions  appear  explicitly  in 
the  integral  equation  (22). 

The  fields  GV  and  u  are,  as  noted,  space-like, whereas  r  a  and 

both  have  time-like  components.  Because  of  this  mix  of  space-like  and  general 

four  dimensional  tensors  in  the  bilinear  form  J^,  the  reduction  of  the  final 

P 

Integral  in  (22)  to  a  surface  integral  over  the  boundaries  of  SI  is  most 
simply  accomplished  by  separating  the  terms  into  pure  space- like  and  time¬ 
like  conponents .  Thus,  writing  out  the  tensors  in  component  form  and 
separating  the  result  into  time  and  space-like  integrals  gives  for  the  final 
integral  in  (22) : 

fh*  [CaTaB  -  “a  " 


(26) 


The  first  of  these  can  be  reduced  by  an  application  of  the  divergence 
theorem  to  the  spatial  coordinates,  so  that: 
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A 


where  the  sum  is  over  the  ordinary  spatial  (subregion)  boundaries 
and  the  integrals  over  time  and  space  are  separated.  We  can  separate  the 
set  of  boundary  surfaces  3v^*^  into  two  groups,  that  is,  into  those 
defined  to  be  the  external  boundaries  of  the  entire  medium  and  those  that 
correspond  to  internal  boundaries.  The  latter,  of  course,  always  are  the 
common  boundaries  between  two  subregions.  Because  the  normals  to  the 

surfaces  are  always  directed  outwardly  from  any  subregion  then  on  an 

interior  boundary  we  always  have  =  -  n^+^ .  Therefore  the  surface 


integrals  arising  from  the  individual  subregion  integrations  on  common 
(internal)  boundaries  can  be  combined,  and  we  obtain: 
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/  >4  [Vkt  ’  Uk^‘] " x  "•[dt°[/5 


Vkt  ~  “k6kt  >n£da, 


Y 


091. 


3v 


(27) 


where  the  double  bracket  notation  in  the  second  integral  over  all  internal 
boundaries  3v^  within  the  medium  denotes  the  change  or  jump  in  the  quantity 
within  the  bracket  across  the  boundary  surface.  That  is. 


H  ■  ^<p)>  -  '><p+1)> 

for  all  p  -  1,2,...,N.  Here  denotes  the  tensor  product  in  the  integrand 
and  J*?(3v^^)  means  the  quantity  evaluated  on  the  internal  boundary 
separating  the  p—  and  p+1—  subregions  approached  (in  the  limit)  from  within 
the  p^  subregion,  while  J^(3v^p+1^)  is  the  same  quantity  evaluated  in  the 
limit  by  approaching  the  boundary  from  the  (p+1)  subregion. 

In  (27)  we  have  suppressed  the  subregion  indexing  (over  p) ,  so  that 
the  sum  of  surface  integrals  over  each  of  the  internal  boundaries  is  replaced 
by  the  convention  that  3Vj  represents  all  the  internal  surfaces..  The  surface 
integration  is  to  be  taken  over  all  the  disjoint  surfaces.  The  same  convention 
of  course,  applies  to  3v„,  which  represents  all  the  external  spatial  surfaces 


of  the  medium. 


The  final  integral  in  (26)  can  be  written 


where  the  spatial  volume  v93v^  excludes  the  internal  boundary  points  and 
implies,,  instead,  limits  from  opposite  sides  of  these  surfaces  approached 
from  within  the  separate  spatial  subregions.  Further,  vG3Vj.  includes  the 
possibility  of  moving  boundaries,  either  external  ox  internal.  Therefore 
v03Vj.  may  be  explicitly  a  function  of  time  t^. 

Following  Archambeau  and  Minster  (1977),  we  observe  that  when 
vGSVj.  is  a  function  of  time  then  the  transport  theorem  provides,  in  effect, 
the  means  of  interchanging  the  partial  differentiation  with  respect  to 
time  and  the  time  dependent  spatial  integration.  Specifically,  for  any 
function  (scalar,  vector  or  tensor  component)  of  the  deformation  or  flow 
in  a  medium,  then  (Archambeau  and  Minster,  1977) 


where  v'  is  any  volume  within  the  medium  and  is  a  component  of  the  boundary 
velocity.  If  U  is  the  particle  velocity  of  the  material  points  then  this 
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Is  the  ordinary  Reynolds  transport  theorem.  The  form  above  is  a  generalization 
of  that  theorem.  Ordinarily,  in  linear  elastic  theory,  is  the  particle 

velocity  and  the  last  term  is  neglected  as  being  of  second  order.  However, 
for  rapid  boundary  movement  —  such  as  for  a  failure  surface  boundary  —  this 
term  cannot  be  neglected. 

If  we  apply  this  theorem  to  the  last  integral  term  in  (26) ,  taking 
F  to  be  the  integrand  in  (26),  and  taking  care  to  partition  the  region 
v03Vj  along  the  internal  boundaries,  we  have: 


(28) 
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Here  we  have  used  Che  fact  that  U.  is  continuous  across  all  boundaries,  by 

^  * 

definition,  and  that  the  normals  from  one  subregion  to  the  next  change  sign 
along  their  common  boundaries. 

The  last  two  integrals  in  (28)  can  be  combined  with  the  integrals  in 
(27).  The  first  integral  in  (28)  is  to  be  treated  as  a  Stieljes  integral 
of  the  form 
t+ 

/df  (t  ) 
m  o 
o 

where  f  (tA)  may  have  step  discontinuities  within  (o,t+)  due  to  discontinuities 

m  o 

along  time-like  surfaces  in  0  (Archaobeau,  1968,  1972;  Minster,  1973; 
Archanbeau  and  Minster,  1977). 

Now,  using  (27)  and  (28)  in  (26)  gives: 

ft+ 

-  u  G\1  d^x°  =  /  dt 
a8  a  aB  I  /  o 

I 


I 


JST,gda 

E 


/  dto  /[❖»]'■*  -  /  *[/  JJdV] 

°  v93v_ 


(29) 


where  we  have  combined  surface  integral  terms  in  (27)  and  (28)  and  used  the 
definition  identities  (in  linearized  form) : 
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Ve 


G 


u  „ 
aBn6 


3G™ 

“  p  3t~  U*  ni 
o 


along  with  the  fact  that 


= 


GMt 


a  aS 


-  u_G 


a  aS 


and 


,4_  T  K  Gm!\l 
JA  ~  p  3tQ  ”  Gk  3tQ  J 


Throughout  we  use  the  linearized  representation  of  the  generalized  space-time 
normal  nQ,  which  is: 

p 


~  ^ni»n2»n3*  “ 


where  n^  is  the  ordinary  spatial  normal  to  a  surface. 
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Consequently,  we  may  now  write  (22)  in  the  form: 


4iru  (x)  -  /  p(x  )  f  (x  )  G^Cxjx  )  d^x° 

V—  I  —o  o—o  a  —  — o 

Jli9ZI 

-  Jdc O  /  J6nBdao  -  /tito/|Jene]<lac 


o  3v, 
.t+ 


o  3v, 


/.[/,„] 

o  v©3\> 


(30) 


Now  we  observe  that  spatial  boundary  conditions  on  the  displacement 

.|- 

field  along  "welded"  internal  boundaries  are  given  in  equation  (12)  as  : 


[VrL^Ivs] 

H.vx " 0 


3V, 


(31) 


Thus  in  (30),  we  have  for  the  integrand,  in  the  integral  over  internal  surfaces: 


M„.  ■  M 


3v, 


Ta3n3  Ua 


M ... 


+The  jump  notation  o  j.  for  a  surface  Z  will  often  be  written  simply  as  o. 
with  the  condition  applying  on  any  of  a  set  of  surfaces. 
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If  then  satisfies  the  same  boundary  conditions  over  the  surfaces 
namely: 


(32) 


then 


itself  is  continuous  at  all  internal  boundaries. 


That  is 


and  the  internal  boundary  integrals  all  vanish.  Thus  (31)  and  (32)  are 
"natural"  boundary  conditions  for  a  region  with  welded  internal  boundaries 
of  discontinuity.  Equations  (32)  are  therefore  "proper"  boundary  conditions 
for  the  Green's  tensor,  in  that  they  are  complementary  to  those  of  (31) 
and  cause  the  internal  surface  integrals  to  vanish. 

We  note,  however,  that  (30)  has  been  obtained  without  direct  use  of 

g 

boundary  conditions  on  either  u^  or  G^.  It  is  therefore  applicable  to  any 
setof  (linear)  boundary  conditions  for  u^  and  any  choice  of  boundary  conditions 

g 

for  G^.  It  is  clear,  however, that  a  choice  of  "proper"  boundary  conditions 

g 

for  G^  will  greatly  simplify  (30)  and  that  if  u^  satisfies  natural  boundary 
conditions  (i.e.  those  naturally  appearing  in  (30))  then  additional  reduction 
of  (30)  is  possible. 

We  will  denote  the  particular  Green's  tensor  that  satisfies  the  boundary 

g 

conditions  (32),  on  internal  boundaries,  by  ^(x;*^).  Thus 
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lv\8 

[  “ylvI  •  0 


(32-a) 


and  Chen  (30)  becomes: 


4ttu^  (x) 


a  -f  o*.»i 


A  o 
d  x 


tP,3  o  1 
J;d  x 
4 

v03v. 


(33) 


where  satisfies  the  boundary  conditions  complementary  to  those  of  (32-a). 

The  external  boundaries  of  the  medium  3v_  may  be  defined  rather  arbitrarily 

t 

in  that  we  may  chose  any  boundary  set  to  be  the  external  boundary  of  the  medium 
in  question  and  this  simply  defines  the  region  over  which  (30)  or  (33)  applies. 

In  particular,  we  normally  chose  one  external  medium  boundary  to  be  any  moving 
boundary  so  that  (30)  or  (33)  applies  to  the  medium  on  one  side  or  the  other 
of  this  boundary.  However,  this  choice  is  not  a  necessary  one  and  if  the 
moving  boundary  is  not  taken  as  an  external  boundary  then  it  is  to  be  treated 
as  an  internal  boundary.  However,  such  boundaries  cannot  generally  be 
considered  as  welded  boundaries  and  therefore  the  second  condition  of  (31) 
does  not  apply.  Instead  the  weaker  condition  insuring  conservation  of  mass 
applies,  that  is,  in  linearized  form: 


t’cVk  »- 0 


However,  Che  firsC  condition  in  (31)  always  applies,  since,  it  states  that 
momentum  is  conserved  across  boundaries.  Thus  if  one  of  the  internal 


boundaries  is  a  moving  boundary,  then  on  this  particular  internal  boundary 
(call  it  3v°)  we  have 


ML  -[•:] 


3v 


o  Ta8nB 


Clearly,  if  we  chose  a  Green's  function  such  that  it  satisfies  the  complementary 
condition 


and  the  "welded  boundary"  condition: 


then 


In  the  situation  in  which  3v°  separates  materials  with  differences  in 
physical  properties  and  when  the  boundary  is  not  geometrically  simple,  then 
it  is  impractical  to  generate  the  Green's  function  satisfying  these  boundary 
conditions.  There  are  several  ways  around  this  difficulty  which  allow  at 
least  good  approximate  results.  However,  we  need  not  consider  them  in 
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detail  here  since  our  ultimate  goal  Is  not  to  treat  the  general  moving 
boundary  problem  but  to  obtain  some  particular  results  that  Involve  fixed 
internal  boundaries  with  a  moving  boundary,  if  present,  as  an  external 
boundary. 

Thus,  In  this  development,  ve  will  take  any  moving  boundary  to  be  an 

external  boundary.  The  treatment  of  the  spatial  surface  Integral  terms, 

in  this  case,  is  formally  the  same  as  when  all  the  boundaries  are  fixed. 

With  this  provision,  we  return  to  (33)  and  note  that  appears  in  the 

p  p 

integral  term  over  3Vg.  Thus  for  a  solution  for  u^  -  rather  than  an 
integral  equation  in  u^  -  it  is  necessary  that  we  know  the  value  of  this 
function  over  the  external  surface.  Since 

Jene  "[HaTaB  -  "c/UK 

and  if  only  the  generalized  tractions  x  n„  are  known  or  specified;  say: 

Clp  P 


We 


*  «=  b 

3ve  a 


(34-a) 


then  an  appropriate  Green's  function  is  such  that 


h>6 


(34-b) 


3v 


E 
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in  which  case  the  second  term  involving  uq  vanishes.  Hence,  should 
satisfy  the  homogeneous  complementary  condition*.  Conversely,  if  u^  is  known. 


so 


uo  1 3Vg  "  Ca 


(35-a) 


Then  the  complementary  homogeneous  condition  for  the  Green's  function  is 


H1 


a  3v, 


In  these  two  cases  we  get,  respectively 


(35-b) 


Jen3  3v 


b  HP 
a  a 


(36) 


and 


JBn6  3u 


E 


CaCnB 


(37) 


In  these  cases,  assuming  can  be  found  satisfying  both  (32-a)  and  either 
(34-b)  or  (35-b)  as  is  appropriate,  then  (33)  can  (normally)  be  used  to  obtain 
a  solution  for  u^(x).  Often  it  is  not  possible  to  find  H^  satisfying  all 
these  conditions  on  all  of  the  boundaries.  In  particular  the  conditions 
(34-b)  and/or  (35-b)  present  considerable  difficulties  on  at  least  parts 
of  the  external  boundary.  Again  there  are  various  approximations  that  can 


-t-Note  that  on  a  free  surface,  that  is  an  external  boundary  in  contact  with  a 

P  I 

vacuum,  then  b  «*  0,  a  ™  1,2, 3, 4.  This  causes  Jeng  to  vanish  when 
,  I  °  3vv 


0  on  such  a  boundary. 


*■ A\ 
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be  used  (e.g.  Archambeau  and  Minster,  1977).  However,  in  some  applications 
both  u^  and  T^Og  can  h®  compatibly  given  as  known  functions  on  parts  or  all 
of  the  external  boundaries.  In  this  case  need  not  satisfy  any  additional 
conditions,  other  than  those  of  (32-a) ,  in  order  to  obtain  a  solution  for  u^. 

This  particular  situation  will  be  treated  in  the  following  sections  since  it 
occurs  for  some  particularly  important  problems. 

The  final  integral  in  (30)  or  (33)  corresponds  to  a  generalized  initial 
value  term  or  contribution  and  has  been  treated  in  this  form  by  Archambeau  (1972) 
and  Archambeau  and  Minster  (1977),  as  well  as  in  other  equivalent  forms  by 
Archambeau  (1964,1968).  In  particular  if  the  spatial  integral  term  In  brackets 
is  continuous  in  (o,t+),  then 

o  v03Vj  903Uj. 

Because  of  the  causal  character  of  the  Green's  functions  used,  the  value  at 
tQ  *•  t+  always  is  zero,  so  that  only  the  value  at  the  initial  time  tQ  =  0 
contributes  and  this  involves  the  initial  values  uq(x£,o)  and  3t  ua(x^,o). 

If  these  initial  values  are  zero,  then  the  integral  term  vanishes. 

However,  in  general  the  integral  can  be  discontinuous  at  a  set  of  times 
throughout  the  interval  (o,t+).  In  this  case  the  time  interval  is  to  be 
partitioned  into  time  segments  in  which  the  integrand  is  continuous,  just 
as  was  done  in  the  case  of  internal  boundaries  in  the  regular  spatial  domain. 


Thus, 


for  time  boundaries  at 


then  we  have: 


(38) 


d  v03v,  P  v93\>. 


(M) 

where  the  range  of  summation  over  the  time  boundary  points  is  to  M,if  t  >  t  , 

/V|\ 

or  I^if  t  <  r  ,  where  L  is  the  index  of  the  discrete  time  point  in  the 
set  closest  to,  but  less  than,  t.  The  bracket  notation  is  the  same  as  that 
previously  used  for  jumps  at  boundaries,  except  here  the  boundaries  are  at 
the  time  points  t^P^ . 

We  observe  that,  by  definition: 


Ik 

v93v 


t(P) 

o 


'  148 


[A 

v93v 


t  = 


t<P>- 


•/ 


,p,3  o 
Jjd  x 


v93v. 


(38-a) 


Since  we  are  considering  the  medium  to  be  bounded  by  external  boundaries  that  may 

move,  then  formally  at  least,  v03v^  may  be  different  in  value  at  t^P^-  e 

and  t^p^+  e.  However,  we  must  define  the  medium  in  which  the  Green's 
o 

function  representation  of  the  field  u^  holds  to  be  that  occupying  the  region 
which  has  not  been  traversed  by  the  moving  boundary.  That  is  the  volume  in 
which  the  integral  representation  of  the  field  u^(x)  applies  must  always 
exclude  the  volume  region  swept  out  by  the  moving  boundary  in  any  time 
increment.  Consequently,  the  change  implied  by  (38)  and  (38-a)  applies  to 
a  volume  integral  over  the  minimum  of  v93Vj  at  t^p^-e  and  t q P^+  c.  Hence, 


Ill 


(38-a)  must  be  written  as: 


where 


(38-b) 


Hence,  the  generalized  initial  value  term  of  (38)  is  given  by: 
t+ 


/-[/'HI  - 17  H 


.3  o 
d  x 


(p) 


(39) 


\  .03vt 

P  1 


Since: 


t(P) 

o 


then  it  is  appropriate  to  require  that  the  Green's  function  have  the  continuity 
properties 


(40) 
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Alternately  we  may  have,  when  impulse  forces  are  applied 


b 


(p) 


-<IML 


(p) 


(44-a) 


and  the  appropriate  time  boundary  condition  for  the  problem  is 


M, 


(P) 


=  0 


(44-b) 


By  far  the  most  important  case  in  applications  is  that  given  In  (43-a) 
and  (43-b).  Specifically,  Archambeau  (1968,1972)  and  Archambeau  and  Minster 
(1977)  consider  failure  in  a  prestressed  medium  and  show  that  (43-a)  and 
(43-b)  can  be  used  to  describe  the  radiation  field  arising  from  the  failure 
process.  In  this  case  (43-a)  applies  and  ff  u^  D  ^  corresponds  to  the 
change  in  the  equilibrium  field  due  to  the  creation  of  a  failure  boundary 
(or  its  incremental  growth)  so  that  in  this  application: 


MU -44 


(P)  3toGk 
o 


where  u*  is  used  to  denote  the  equilibrium  field  explicitly. 

Equation  (39)  expresses  the  contribution  to  the  field  u^ (x)  arising 

from  the  discontinuous  behavior  of  the  spatial  volume  integral  of 

|  (  .  )M  4 

at  a  discrete  set  of  times  (t'”^  f  .  We  may  generalize  this  result  to 


include  a  representation  in  which  the  set  <  t 


(P) 


M 


becomes  partly  or  totally 


a  continuous  distribution.  The  result  is  a  straightforward  generalization 
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of  (39)  which  has  been  shown  by  Archambeau  (1968,  1972)  and  Archambeau 
and  Minster  (1977),  in  the  somewhat  more  restricted  context  of  a  growing 


failure  surface  problem,  to  be  of  the  form: 
t+ 


min(t,t^) 
o 


o  v03v  o  v03v 


(45) 


(M) 

where  t  is  the  upper  limit  point  of  the  continuous  and/or  discrete  sec 
o 

of  discontinuity  points  and  6JP  represents  the  incremental  change  or 
variation  of  Jp  over  a  time  increment  6tQ.  In  (45)  all  quantities  are  to 
be  measured  relative  to  the  initial  state.  (See  Archambeau  and  Minster,  1977 
for  other  choices  for  the  reference  state.)  Here  6JP  is  equivalent  to 
(I  U  in  (39) ,  while  the  time  integral  over  t  replaces  the  sum  over  the 
index  p.  Formally,  the  integral  is  the  limit  of  the  summation  as  the  time 
spacing,  6 1, between  successive  t^P^  is  allowed  to  become  infinitesimal. 

To  explicitly  display  the  discrete  discontinuous  case,  it  is  only 
necessary  to  note  that  discrete  discontinuities  in  at  times  t^P^  correspond 
to  step  function  discontinuities  of  magnitude  [[  J  ,  at  tQ  *»  t^P^.  Then 


with  6(t  -t  ^)  the  Dirac  delta  function.  It  is  clear  that  use  of  this  in 
o  o 

(45)  gives  the  discrete  result  (39).  Hence  (45)  is  a  general  form  of  the 

result  Including  both  discrete  and  continuous  distributions  of  time  disasntinuities. 


For  the  cases  of  greatest  interest,  namely  when  3^  u^  is  continuous 

o 

and  G®,  and  its  first  derivative  in  tfl  are  continuous,  then  the  integrand 
in  (45)  has'the  explicit  form  (e.g.  Archambeau  and  Minster,  1977) 


(5J4/Se„)-[P  \>]v” 


(46) 


where  u^(x^.t^)  is  the  equilibrium  field  value  at  the  source  time  t^.  The 

* 

factor  6u^/6to  has  been  written  as  a  partial  derivative  with  respect  to 
t  since  the  meaning  is  the  same  as  the  variation  ratio.  The  general  case 
is 


M  - •> (5v5t„)stoG”  -  < 


where  vk  =  J  V 
o 

To  collect  and  summarize  the  most  useful  results  of  this  section,  we 
have  shown  that  in  general,  the  integral  Green’s  tensor  representation  of 
u,(x)  is: 


4iru^(x) 


!x)  =  /  dto  j  pf 

°  v©3v. 


pf  c“dV 

a  a 


-f  dto  /  W*o  -/  dto/[Js’'6]h 

o  3v_  o-  3v_ 


(47) 
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M 

where  the  result  (45)  has  been  used  In  (30)  and  t  =  min(t,to).  The  Green’s 

tensor  used  to  obtain  (47)  is  the  causal  solution  of  L  Gp  =>  Ap  satisfying 

oy  Y  a 

g 

the  condition  that  Gp  and  its  first  time  derivative  with  respect  to  t  be 

Y  o 

continuous  for  t  <  t.  When  the  medium  internal  boundaries  are  fixed  and 
o 

welded  and  any  moving  boundaries  are  external  boundaries,  then  u^(x)  satisfies 
the  boundary  condition 


[V,]3Vi  5  [We] 

M  “  ° 


3v 


(48) 


on  all  internal  boundaries.  Then, 


IM3V  ■  UL 


(48-a) 


Further  for  spontaneous  processes  giving  rise  to  generalized  initial  values, 
such  as  failure  processes  involving  moving  (or  growing)  boundaries  in  an 
initially  stressed  medium,  then 


(5J«/6t)  -  _» 8to\]  \Gl 


with  u^  the  equilibrium  field  for  the  medium  within  the  external  boundaries. 

*  * 

Here  3fc  u^  is  the  "source  time"  derivative  of  u^  and  is  the  source  term  for 

°  * 
the  radiation  field  arising  from  the  process.  In  principle  u^(x^,tQ)  can  be  obtained 
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Independently  from  (47)  (e.g.  Archambeau,  1968)  so  that  the  integrand  of  the 

final  term  in  (47)  can  be  considered  as  specified.  Likewise  the  applied 

force  field  f  (x  )  can  be  considered  as  a  known  field.  Thus  these  contributions 
o—o 

to  can  be  obtained  directly,  assuming  that  G^, satisfying  the  conditions 

Just  stated,is  obtained.  In  view  of  (48-a)  however,  the  surface  integral  over 

3v_,  the  internal  boundaries,  contains  u  so  that  (47),  as  it  stands,  is 
I  a 

an  integral  equation  for  u  .  Further  the  surface  integral  over  3v  ,  the 

o  £ 

external  boundaries,  may  contain  u^  terms  as  unknowns. 

If  we  use  the  particular  Green's  tensor  H^(x;x^  which  satisfies  the 
internal  boundary  conditions  complementary  to  those  for  uq  in  (48),  namely 

14  • 0 

then  with  HW  used  for  G*1  in  (48-a)  we  have 
a  a 

w9V- 

and  the  integral  over  the  internal  boundaries  in  (47)  vanishes. 

Further,  if  the  external  boundary  of  the  medium  is  partly  or  wholly  a 
free  surface,  then  the  total  external  boundary  may  be  denoted  as 


3ve 


+  3v 


(1) 

E 


(50) 
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with  denoting  a  free  surface  and  the  remaining  part  of  the 

£  £ 

external  boundary.  On  3\>g^  the  boundary  condition  for  u^  is: 


3v 


(0) 


■  0 


(51) 


The  complementary  condition  for  the  Green's  tensor  is: 


fTc  „ph 

„  Tl 

|[8cnrHY 

HoB"B 

L  JJ. 

3v 


=  0 
(0) 


(52) 


In  this  case,  using  H ^  satisfying  (49)  and  (52)  with  u^  satisfying  (48) 
and  (51) 


Thus,  with  the  Green’s  tensor  H^,  the  integral  relation  for  u^(x)  is: 


t+ 

4iru  (x)  -  f  dt  /  pf  HV  d^x° 
li  —  /  o  /  a  a 

o  vOSVj 

“A/  JeVao 

/ 1  \ 


(53) 


3v 


(1) 
E 


/*•/[’  v«‘KH«d: 

o  v03v. 
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where 

[  Jen8  =  Va3  "  UaHae]  n3 

If  values  of  the  fields  u  or  t  0n0  are  known  or  specified  on 

a  aS  3  E 

then  complementary  homogeneous  conditions  on  or  reduces  on 

to  a  known  function,  not  involving  the  unknown  u^  or  its  derivatives,  and 

(53)  is  a  solution  for  u^(x).  That  is, in  this  case,  in  combination  with 

its  space  and  time  derivative, act  as  transfer  functions  for  the  field 

specified  on  space  and  time  boundaries, such  that  the  field  is  propagated  to 

other  spatial  points  at  later  times. 

Alternately,  if  compatible  values  of  both  uq  and  T-flg  are  known  or 

can  be  specified  on  ,  then  need  not  satisfy  any  additional  condition 

in  order  that  (53)  provide  a  solution  for  u  (x).  Th-U.  is  J„n„,  in  this 

y  p  p 

case,  is  a  known  function  on  when  the  Green's  tensor,  H^1,  satisfying 

(only)  (49)  and  (52)  is  used. 

In  the  following  sections  this  latter  case  will  be  demonstrated  to 
occur  in  several  important  applications.  In  order  to  use  (53),  however,  must 
be  specified.  The  generation  of  for  media,  of  the  greatest  interest  in 
elastodynamic  theory  applications  to  Geophysics  are  layered  elastic  spheres 
and  half  spaces.  The  form  of  in  such  media  will  therefore  be  considered 
in  following  sections. 
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(4)  Green’s  Tensor  Integral  Representations  in  the  Frequency  Domain 

For  this  development  to  be  pursued,  it  is  necessary  to  transform  all 
the  equations  and  fields  from  the  time  domain  to  the  frequency  domain.  This 
was  the  approach  used  by  Archambeau  (1968)  from  the  onset  in  treating  the 
particular  problem  of  a  growing  rupture  zone  in  a  stressed  medium,  this 
being  a  particular  case  of  the  general  problem  treated  here.  We  can  however 
use  the  time  domain  integral  representation  for  uq  given  in  the  previous 
section,  to  obtain  the  equivalent  frequency  domain  result. 

In  particular,  we  define  the  Fourier  transform  operation  with  respect 
to  time  t  by 

“«  2  4  e_1“dC  <54> 

—00 

Here  denotes  the  operator  and  u^  the  transformed  function.  The  inverse 
operator  ^  ^  is  defined  by 


where  u>  ■  2irf  is  the  transform  variable  and  here  corresponds  to  angular 
frequency  with  f  the  frequency. 

The  convolution  of  two  (tensor)  functions  is  then  given  by 

i  1  { vJ  =  /  u  (t  )  v  (t-t  )dt 
u  v  Q  p;  J  a  o  g  o  o 

— OO 
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and  If  u  is  such  that  u  (t  )  =  0  for  t  <0  while  v  (t  )  =  0  for 
a  a  o  o  —  a  o 

t  ^  t  =  t  +  e,  then  this  becomes: 
o  * 

t+ 

ll  )  I  ”  /  u  (t  )vc(t-t  )dt 

a<i>  (  a  3  /  J  a  o  $  o  o 

o 

In  this  case 

r t+ 

it  Ve}  s  Vb  ’  (t{J  “«(to)vB(t-to)'J':o } 

o 

Convolution  relations  of  this  type  can  be  seen  to  appear  in  the  Green’s 
integral  relations  for  the  elastic  field  u^Cx),  for  example  in  (53).  Thus 
since  Vg  in  the  above  relations  plays  the  role  of  a  transfer  function,  then 
it  is  clear  that  the  various  integral  kernels  involving  Green's  functions 
play  the  same  role.  This  will  become  even  more  evident  in  the  following 
development . 

Focusing  our  development  on  the  most  useful  of  the  Green's  integral 

relations,  we  consider  the  Fourier  transform  of  uq(x)  in  (53)  with  respect 

to  the  (observers)  time  variable  t.  In  order  to  carry  out  this  operation 

in  the  most  routine  way,  we  first  observe  that  the  two  point  Green’s  tensor 
p  lc  1c 

Hq(x  ,t;xo,tQ)  is  causal,  so  that 

HW(x;x  )  =  0,  for  t  >  t 
a  — o  '  o 

This  allows  us  to  extend  the  integration  range  in  the  first  and  second  integral 
terms  from  (o,t+)  to(o,®)  without  chanlng  the  value  of  these  integrals. 
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since  H*1  vanishes  identically  when  t  >  t+. 
a  o  ~ 


Next  we  note  that: 


u  (t  )=  0  for  t  <  0 
a  o  o 


by  supposition  (i.e.  the  field  effects  start  at  t  =  0).  This  allows  us  to 

integrate  over  the  "source  time"  range  (-“,+«»)  in  the  first  two  integrals 

since  these  integrals  vanish  in  the  range  (-“,o)  because  of  the  latter 

condition.  Finally,  for  the  last  integral  Involving  the  "generalized  initial 

values",  we  can  take  the  time  changes  in  the  equilibrium  field  to  be  zero 

M 

in  the  range  outside  the  interval  (o,tQ)  since  such  changes  only  occur  in 
this  interval  by  our  earlier  definitions.  Thus 


u*(t  )  »  0;  t  <  0,  t  >  tM 
to  o  o  ’  o  o 

o 


This, in  combination  with  =  0  for  t  >  t,  shows  that  the  integration  range 
for  the  final  integral  can  be  extended  to  (-",+*“)  without  changing  the  value 
of  the  integral  as  well.  Thus  (53)  can  be  written  as 


(55) 
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Finally  we  observe  that  the  Green's  tensors  are  functions  which  depend 
only  on  the  time  difference  t-to.  This  follows  from  the  fact  that  they 
are  required  to  be  causal  and  because  of  the  self-adjointness  of  the 


generating  differential  operator  L^.  Thus  we  have  that: 

Hy  =  HW(xk,xk,t-t  ). 
a  a  o  o 

Now  operating  on  equation  (55)  with  where  we  note  that  this  operator 

commutes  with  all  the  integral  operators  on  the  right  side  of  (55),  we  get: 
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If  the  medium  boundary  3v^  is  not  stationary,  then  vQSv^  and  will  be 

functions  of  t  ,  so  that  the  space  integrals  in  (56)  must  be  evaluated 
o 

first  and  then  the  integral  over  tfl  may  be  obtained.  Thus  the  time  integral 

over  t  which  is  itself  in  the  form  of  a  Fourier  transform  operator,  does 
o 

not  commute  with  spatial  integral  when  there  is  boundary  movement.  (Of 
course  approximations  can  be  obtained  by  neglecting  the  volume  or  surface 
area  changes  and  treating  the  limits  of  the  spatial  integrals  as  fixed  in 
time,  at  least  to  first  order.) 

If  all  the  boundaries  of  the  problem  are  fixed  (in  the  sense  of  the 
usual  elastic  theory  approximation  in  which  boundary  movement  with  the 
particles  is  considered  as  a  fixed  boundary) ,  then  the  integration  with 
respect  to  t  and  the  space  integrals  can  be  interchanged.  We  have  then: 


4xu  ”  /  pf  H^d^x 
y  J  M  a  a  c 


v93v. 


-  /  [  fiWT  -  U  1 

J  L  a  aB  a  aP  J 


V*o 


3v, 


(1) 

E 


(57) 


-  f  pa*  H^d^x 
j  a  a  c 


v03v. 


Here  the  boundaries  are  to  be  considered  as  fixed  after  the  time  t  *  0. 

o 

However,  this  representation  still  admits  the  case  of  an  instantaneously 

*  k 

created  boundary  at  tQ  •  0.  In  this  case  u^x^.t^)  is  a  step  function  at 

t  ■  0  whose  transform  u*  is  proportional  to  l/u>.  This  problem  was  treated 
o  o 
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by  Archambeau  (1972)  using  a  representation  form  of  essentially  the  type 
(57),  with  only  the  last  term  retained  as  significant. 

It  is  also  worth  pointing  out  that  since  the  boundaries  are  assumed 
fixed  here,  then  =  0  in  (57).  In  this  case  all  the  tensors  are  purely 
space- like;  for  example  =  (n^,^ ,n^,0) .  This  means  that  all  the  Creek 
indices  in  (57)  may  be  replaced  by  three  term  indices,  denoted  by  the  latin 
forms,  i,  k,..,  etc. 

The  general  structure  of  (57),  when  viewed  in  terms  of  the  convolution 
definitions  previously  stated,  is  that  of  a  sum  of  convolution  operations. 

In  particular,  all  the  terms  in  (57)  involving  forces,  that  is  pf^,  Tagr,g 

2-*  ~n 

-pw  u^,  have  as  a  (time-like)  transfer  function.  On  the  other  hand, 

the  term  involving  the  boundary  displacement  u  on  ,  has  HP n„  as  a 

a  b  op  p 

transfer  function.  Moreover,  the  Green’s  function  will  be  dependent  on  the 

lc  lc  |( 

spatial  coordinates  x  and  xq  through  the  absolute  difference  x  -x^.  That 

is,  its.  functional  dependence  can  be  expressed  as 


HV  =  HP  (|xk-xk|,  t-t  ) 
o  o  '  o'  o 


This  follows  from  the  causal  condition  imposed  on  H^,  coupled  with  the 

3  6 

character  of  the  generating  equation  l^ir  =  These  of  course  are  the 

same  conditions  leading  to  the  time  difference  dependence  used  earlier  and 

lead  to  the  usual  statement  of  reciprocity  for  the  Green's  function. 

Therefore  with  HP  of  the  functional  form  described,  it  is  clear  that 
a 

the  spatial  integrations  in  (57)  are  also  in  the  form  of  (three  dimensional 

space)  convolutions.  Therefore  HP  and  HVon0  are  seen,  from  (57),  to  have 

a  ap  p 

the  character  of  space-time  transfer  functions  or  propagators. 
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The  form  of  (57)  is  considerably  simpler  than  that  of  the  more 
general  result  (56),  which  admits  of  the  treatment  of  moving  boundaries. 
However,  the  structure  of  the  two  integral  representations  is  essentially 
the  same.  In  particular,  H^  and  are  the  propagators  for  the  fields. 

Since  is  related  to  within  the  integrands  in  (56)  and  (57)  by: 


3Hk 


HWQ(x;x  )=*C  ,(x)  . 

aS  —  —o  aSyfi  -o  3x6 

o 


Then  the  time  transforms  with  respect  to  t  =  x  are  related  by 


C  _  .  (x  ) 
agy6  -t> 


afi11 

3x6 

o 


(58) 


Therefore  knowledge  of  H^,  the  transform  with  respect  to  the  receiver  time 
coordinate.  Is  sufficient  to  specify  the  tensor  appearing  es  a  transfer 
function  in  both  (56)  and  (57). 

Thus ,  it  has  now  been  shown  that  knowledge  of  is  all  that  is  required 
in  order  to  use  (56)  and  (57)  to  solve  a  large  class  of  elastodynamic  problems, 
including  problems  involving  moving  boundaries,  as  well  as  classical  space- 
time  boundary  value  problems.  Therefore  solutions  to  a  number  of  important 
elastodynamic  problems  require  the  solution  of  the  receiver  time  transformed 
equation  corresponding  to: 


L  Hp(x;x  ) 
ay  y  ~  -o 


&B(x;x  ) 
&  —  — o 


(59) 
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with  boundary  conditions,  from  equations  (49)  and  (52), 


(60) 


By  our  previous  definitions  of  the  boundaries  8v  and  ,  These  are  all 

X  tj 

fixed  boundaries  while  only  may  be  a  moving  boundary.  Thus  all  the 

tensors  in  (60)  are  space  like  (i.e.  ,  =  0  on  these  boundaries)  and  so 

these  conditions  may  be  written  in  terms  of  space-like  components  above. 

We  have: 


(60-a) 


Taking  the  transform  of  these  equations  for  IT  gives: 
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vs + »“2  s;  -  -  k  s<*k  -«o> 


where  is  a  function  of  the  ordinary  space  coordinates  alone ,  and  where 


l  «.  1 _  (Q  - _ ) 

Ik  3X,  v  £jkm  dx  ’ 

i  J  m 


is  the  regular  elastic  operator.  The  boundary  conditions  become: 


IM 


M 


where 


H“  - 


k  "  Ckttj  3Xj 


The  boundary  conditions  may  be  written  in  a  contact  matrix  form, 
analogous  to  that  expressed  earlier  in  the  equation  (14),  as: 


[CM 


(62-a) 


i 
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[8lk^]3v«» 

E 

where  B^  is  the  "space  boundary  operator" 
B*k  “  “j  [°£jki  3^] 
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